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analysis.  In  part  III  we  combine  the  ideas  presented  in  the  first  two  parts,  in  our  study  of 
accuracy  stability  and  convergence  of  the  spectral  Fourier  approximation  to  time-dependent 
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1.  SPECTRAL  APPROXIMATION  THEORY 


1.1.  The  Periodic  Problem  -  Fourier  Approximation 
Consider  the  first  order  Sturm-Liouville  (SL)  problem 

(1.1.1)  —<j>  =  \<j>(x),  0  <  x  <  27 r, 

augmented  by  periodic  boundary  conditions 

(1.1.2)  m  =  tw¬ 
it  has  an  infinite  sequence  of  eigenvalues,  A*,  =  ik,  with  the  corresponding  eigenfunctions 
tk(x)  —  etk=-  Thus,  (Afc  =  ik,  <J) *.  =  etkx )  are  the  eigenpairs  of  the  differentiation  operator 
D=j-  in  L2[ 0, 2 7 r),  and  they  form  a  complete  system  in  this  space.  This  system  is  complete 
in  the  following  sense.  Let  L2[0,2tt)  be  induced  with  the  usual  Euclidean  inner  product 

/•2ir  _ 

(1.1.3)  (w1(x),w2(x))  =  /  w1(x)w2(x)dx. 

Jo 

Note  that  <f>k{x)  are  orthogonal  with  respect  to  this  inner  product,  for 

<1X4>  (e’“'e'>)  =  {  lle-llU*  \  Vk. 

Let  tu(a;)eL2[0,27r)  be  associated  with  its  spectral  representation  in  this  system,  i.e.,  the 
Fourier  expansion 

(1.1.5)  w(x)  ~  f)  w(k)Mx)>  Hk)  = 
or  equivalently, 

00  1  /•2tt 

(1.1.6)  w(x)  ~  w(k)elkx,  w(k)=—  w(()e~l'r\ 

fc=-oo  ^ 

The  truncated  Fourier  expansion 

N 

(1.1.7)  SNw=  £  w(k)eikx, 

k=-N 


1 


denotes  the  spectral-Fourier  projection  of  w(x)  into  7r^-the  space  of  trigonometric  polyno¬ 
mials  of  degree  <  N: 

N 

SNw  =  w(0)  +  £[t2>(£)e,7:i  +  w(-k)e~ikx] 

h=l 

N 

(1-1.8)  =  w(0)  -f-  +  rh(— A:)]  cos  kx  -f  —  w(— fc)]  sin  kx 

k=  1 

N 

=  ^2  '^k  cos  kx  -f  bk  sin  kx ; 
k=  0 

here  ak  and  bk  are  the  usual  Fourier  coefficients  given  by 

1  r2* 

a,k  =  w(k)  +  w(- k)  =  —  /  u>(£)  cos  k£d£, 

it  Jo 

(I-1-®) 

Ik  —  i[w(k)  —  w{—k )]  =  —  [  w(£)s\nk£d£. 

7T  Jo 

Since  w  —  S^w  is  orthogonal  to  the  Tr^r-space: 

(1.1.10)  (w  -  SNW,e,,cx)  =  2irw(k) -2irw(k)  =  0,  |A;|  <  N, 

it  follows  that  for  any  pn^n  we  have  (see  Figure  1) 

(1.1.11)  ||«>  -  Pn\\2  =  ||w  -  5'iv^ll2  +  ||Sjvu>  -  Pat||2. 

Hence,  S^w  solves  the  least-squares  problem 


w 


(1.1.12)  ||u»  -  SNw ||  =  MinPN£7rjv||iv  -  pw|| 
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i.e.,  S?jw  is  the  best  least-squares  approximation  to  w.  Moreover,  (1.1.11)  with  p^  =  0  yields 


(1.1.13) 


HStfU/H2  =  ||u/||2  —  \\w  —  i?^to||2  <  ||u/||2 


and  by  letting  N  — >  co  we  arrive  at  Bessel’s  inequality 


(1.1.14) 


2*  2  |»(fc)|2  =  Y,  l*MI2IMI2  <  IMI2- 

jfc=  — OO  Jc=  —  oo 


Remark:  An  immediate  consequence  of  (1.1.14)  is  the  Riemann-Lebesgue  lemma,  asserting 
that 

w(k)  =  —  /  w{()e~lk^d( — *  0,  for  any  weL2{ 0,27r). 

Z7T  JO  k-+oo 

The  system  {<f>k  =  e*fcx}  is  complete  in  the  sense  that  for  any  w(x)eL2[0, 2n)  we  have 
Parseval  equality: 


(1.1.15) 


2*  Y  l“MI2  =  E  K*)I2II«2  =  IMI2, 


which  in  view  of  (1.1.13),  is  the  same  as 


(1.1.16) 


Jim  || to  =  ^  w(k)e,kx  —  tu(x)||  =  0. 

-N 


The  last  equality  establishes  the  L 2  convergence  of  the  spectral-Fourier  projection,  Snw(x), 
to  to(x),  whose  difference  can  be  upper  bounded  by  the  following 

Error  Estimate: 

iiu.  -  =  mi2  -  ns*«,|i!  =  E  \mnh\\2=i*  y  i*wi2- 

|*|>JV  \k\>N 

We  observe  that  the  RHS  tends  to  zero  as  a  tail  of  a  converging  sequence,  i.e., 

,2*  N 

(1.1.17)  /  j-tx;(cc)  —  ^  w(k)eikx\2dx  =  2tt  ^  |u>(A;)|2 — >  0. 

J°  k=-N  |fc|>W  N-*°° 

The  last  equality  tells  us  that  the  convergence  rate  depends  on  how  fast  the  Fourier  coeffi¬ 
cients,  w(k),  decay  to  zero,  and  we  shall  quantify  this  in  a  more  precise  way  below. 

Remark  (1.1.17)  yields  uniform  a.e.,  convergence  for  subsequences;  in  fact  one  can  show 


(1.1.18) 


a.e.  Jim  |to(x)  —  5Arpu;(a:)|  =  0,  inf  >  1. 
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In  fact,  w(x)  =  a.e.  lim^oo  Snw(x)  for  all  weL2[0, 2i r],  but  a.e.  convergence  may  fail  if  to(-) 
is  only  ^[O^Trj-integrable. 

Yet,  if  we  agree  to  assume  sufficient  smoothness,  we  find  the  convergence  of  spectral- 
Fourier  projection  to  be  very  rapid,  both  in  the  L 2  and  the  pointwise  sense.  To  this  we 
proceed  as  follows. 

Define  the  Sobolev  space  Hs[0,2tt)  consisting  of  27r-periodic  functions  which  their  first 
s-derivatives  are  X2-integrable;  set  the  inner  product 


(1.1.19) 


*  [2*  _ 

(wi, w2)h‘  =/,  I  Dpwi(x)Dpw2(x)dx. 

p=cr0 


The  essential  ingredient  here  is  that  the  system  {e’fcx}  -  which  was  already  shown  to  be 
complete  in  T2[0,2tt)  =  H°[ 0, 27r),  is  also  a  complete  system  in  H’[0,2n )  for  any  s  >  0.  For 
orthogonality  we  have 


(1.1.20) 


0  j 

2-k  E  k2p  j  =  k. 

p=0 


The  Fourier  expansion  now  reads 

(1.1.21) 


Vi(x)  ~  Y.  w,(k)c' 


where  the  Fourier  coefficients,  w,(k),  are  given  by 


(1.1.22) 


_  (M(x),e-fa)„. 

>  ~  (eikx,  eikx)nt  ' 


We  integrate  by  parts  and  use  periodicity  to  obtain 


(w(x),  e'kx)ffi  =  E  f  Dpw(x)Dpelkxdx  = 
p=cr° 

=  E(-!)P  /  w(x)D2p exkx dx 
p= o  Jo 

3  f2ir 

=  E  {-m-W 

p=  0  Jo 

and  together  with  (1.1.20)  we  recover  the  usual  Fourier  expansion  we  had  before,  namely 


(1.1.23) 


w3(k)  =  w(k )  =  f  w(()e~lHd(. 
Ztc  Jf= 0 
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The  completion  of  {etfcr}  in  H3[0,2n)  gives  us  the  Parseval  equality,  compare  (1.1.15) 


(1.1.24) 


Since 


(1.1.25) 


||w -&**.=  £  wmle*! H-=  £  ll*MI2  •  2*  ]T>2,>]  > 

|Jt|>iV  \k\>N  P= o 

J2N2p-2t-  £  |tl;(A:)|2  =  ^^2p-||u,-5wu;||2. 

p=0  |A:|>W  p=0 


Const1(l  +  N2y/2  <  f  £  N2pJ  <  Const2(l  +  IV2)?, 


we  conclude  from  (1.1.24),  that  for  any  toeiiP[0,27r)  we  have 


(1.1.26) 


\w  —  5/srtv||  <  Const,  ■  — ,  u>eiV*[0,27r). 


Note  that  Const,  =  Consti  ■  j|u;  —  — *  0-  This  kind  of  estimate  is  usually  refered 

N-*  OO 

to  by  saying  that  the  Fourier  expansion  has  spectral  accuracy,  i.e.,  the  error  tends  to  zero 
faster  than  any  fixed  power  of  N,  and  is  restricted  only  by  the  global  smoothness  of  w(x). 

We  note  that  as  before,  this  kind  of  behavior  is  linked  directly  to  the  spectral  decay  of 
the  Fourier  coefficients  in  this  case.  Indeed,  by  Cauchy-Schwartz  inequality 

KJfe)|  =  h&,(Jfc)|  <  < - i - rlMI*. 

v  ;  Wl  "  i \e'kx\\2H-  -  (27T  0*2p)’ 


(1.1.27) 


<  Const  • 


(HI#' 


In  fact  more  is  true.  By  Parseval  equality 


IN*.  =  £  £  £** 

/c=  —  oo  k=z  —  oo  \p=0  J 

and  hence  by  the  Reimann-Lebesgue  lemma,  the  product  (l-J-|/t|2)?|u>(&)|  is  not  only  bounded 
(as  asserted  in  (1.1.27),  but  in  fact  it  tends  to  zero, 

(l  +  IA^KAOl— >°- 

k-+oo 

Thus,  i u(k)  tends  to  zero  faster  than  jk|~s  for  all  w[x)eH3 .  This  yields  spectral  convergence, 
for 

II™  -  $HI2  =  2?r  E  K*0I2  ^  Const-  E  (TT 7172-77  ^  Const- 
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i.e.,  we  get  slightly  less  than  (1.1.26), 

||n;  —  <Sjvn>)|  <  Const. - - — *  0  s  >  1. 

Ns~2  N—oo 

Moreover,  there  is  a  rapid  convergence  for  derivatives  as  well.  Indeed,  if  w(x)eH3[ 0, 27 r) 
then  for  0  <  a  <  s  we  have 

||u>  -  Snw\\2h.  =  •£(2','£,k*)\w(k)\’ 

|fc|>AT  p=  0 

<  Const.  (1  +  |A;|2)<r |t£>(A:)|2 

<  Const.  V  < 

£h  (1  + 

£  Con,t-1£  (1+^)-  w )!  = 


<  Const. 


II w  -  SaHIh* 

N*-”) 


Hence 


(1.1.28) 


w  —  ^HIk0,  <  Const,  •  — — • -,  a  <  s,  weH3[0,2%) 


with  Const,  ~  \\w  —  SVrojltf. — *  0.  Thus,  for  each  derivative  we  “lose”  one  order  in  the 

N-*oo 

convergence  rate. 

As  a  corollary  we  also  get  uniform  convergence  of  S^w(x)  for  i?1[0,27r)-functions  w(x ), 
with  the  help  of  Sobolev-type  estimate 

(1-1.29)  Maxo<x<2*Kx)|  <  Const.||x>||iyi . 

(Proof:  Write  v(x )  =  v(xo )  +  J^v'(x)dx  with  u(so)  =  i*fo*v(x)dxi  and  use  Cauchy- 
Schwartz  to  upper  bound  the  two  integrals  on  the  right.) 

Utilizing  (1.1.20)  with  v(x)  =  w(x)  —  S^w(x)  we  find 

Maxo<x<2*Ks)  -  57/ru(®)|  <  Const. ||u.  - 

(1.1.30) 

<  Const,  — — — — *  0,  uieif^O,  2tc). 

1  N-*  oo 

Corollary:  Assume  w(x)eH1[ 0,2tt).  Then 


(1.1.31) 


w(x)  =  ^2  w(k)e,kx. 


6 


JLI^hItrli°n’  "  n°te  ‘hat  ^  ““  *“*<*».  can  be 


(1.1.32a) 


where 


(1.1.326) 


fc=-AT 

=  L0Dn(x  -  (MOd( 


dn(x  -  ( )  =  L  y  ei^-si  =  i_ sil1  (w  +  5)  (*  -  f) 

w  2-  sin  (S)  • 

Thus,  the  spectral  projection  is  given  by  a  convolution  with  the  so-caUed  Dirichlet  kernei, 

(1.L33)  D,.  =  J.fi^+|)x 

2-n-  '  'r 

Now  (1.1.23)  reads 


smf 


(1'U4)  Wl)  -  *  **  *  ■  ~  Const,  ~ 
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1.2.  The  Pseudospectral  (Collocation)  Approximation 


We  have  seen  that  given  the  “moments” 

w(k)  =  r  w{i)e~ikid(,  -N  <  k<  N, 

2?r  J£=  0 

we  can  recover  smooth  functions  w(x)  within  spectral  accuracy.  Now,  suppose  we  are  given 
discrete  data  of  u/(x):  specifically,  assume  to(x)  is  known  at  equidistant  collocation  points 

(1.2.1)  wl/  =  w(xu),  xv  =  r-\-vh,  v  =  0, 1,-  •  •  ,2N. 


Without  loss  of  generality  we  can  assume  that  r-which  measures  a  fixed  shift  from  the  origin, 
satisfies 

(1.2.2)  0  <r<h  =  ------ 

v  '  2N  + 1 

Given  the  equidistant  values  wu ,  we  can  approximate  the  above  “moments,”  w(k),  by  the 
trapezoidal  rule  2 


(1.2.3) 


?W+1 


™(*0  =  2^  £  "w»e 


"...  _ 


i/=0 


2iV  + 


2N 

T£ 

1  i/=0 


wve 


—ikxu 


Using  w(k )  instead  of  w(k )  in  (1.1.7),  we  consider  now  the  pseudospectral  approximation 


N 


ipNw  =  w(k)e'kx. 
k=-N 


(1.2.4) 

The  error,  tn(x)  —  ip^w(x),  consists  of  two  parts: 


w(x)  -  ipNw(x)  =  w(k)e,kx  -fi  [u>(A:)  -  w(k))eikx. 

|A=|>AT  |*|</V 

The  first  contribution  on  the  right  is  the  truncation  error 


(1.2.5)  Tnw(x)  =  (I  -  SN)w(x)  =  J2  w{k)eikx. 

\k\>N 

We  have  seen  that  it  is  spectrally  small  provided  w(x)  is  sufficiently  smooth.  The  second 
contribution  on  the  right  is  the  aliasing  error 

(1.2.6)  A^w(x)  =  [tw(fc)  —  w(k)]elkx. 

\k\<N 

This  is  pure  discretization  error;  to  estimate  its  size  we  need  the 

2£/  and  indicate  summation  with  |  of  the  first  and  respectively,  the  first  and  the  last  last  terms. 
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Poisson’s  Summation  Formula  (Aliasing) 


Assume  to(a;)ei/’1[0, 2w).  Then  we  have 

CO 

(1.2.7)  w(k)=  Y,  eip^N+1)rw(k  +  p(2N+l)). 


p—  —  co 


Proof:  For  w(x)eH1[ 0,2tt)  we  insert  its  Fourier  expansion 


(1.2.8) 


1 


2  N 


—ikx  v  _ 


i  2W 

2JV  +  1~0 


E 


j=— oo 


—  t/cav 


Since  u/(a:)  is  assumed  to  be  in  Hl,  the  summation  on  the  right  is  absolutely  convergent 


(1.2.9)  E  W)\  ^  (B1  +  ;2Mi)!2  '  E  YJpj  <  Const. Him, 


(1.2.10) 


1  co  2N 

*(*)  =  WTT  E 

<i-,v  ■*  1  j'=-oo  i^=0 


and  hence  we  can  interchange  the  order  of  summation 

2  AT 

+ 

Straightforward  calculation  yields 

2N 

L 

2  N 


i  2N  i  2N 

1 _ V'  Ai-^lr+uh)  _  .  1  .  JU-k)v*T%T  — 

^  +  l^o  2AT4-  1 


(1.2.11) 


=  e i{j~k)r  ■ 


1 


(  e«j-k)l*JlT+ll  _  i 

2Ar+x  =0  j  —  k  ^  0[mod  2Ar  +  1] 


2iV  +  1 


e‘k'-fc)5WTr  _  i 


[  2N  +  1, 

and  we  end  up  with  the  asserted  equality 

2  N 


j  —  k  =  p  •  ( 21V  -f  1). 


w(k) 


CO  1  JiV  CO 

=  E  *«)  •  ^rrr  E  =  E  *(«=  +  p(2At  + 1))  • 

j—  —  oo  1  i/=0  p=-oo 


We  note  that  once  w(x)  is  assumed  to  be  smooth,  it  is  completely  determined  (pointwise, 
that  is)  by  its  Fourier  coefficients  w(k)\  so  are  its  equidistant  values  wv  =  w(xl/)  and  so  are 
its  discrete  Fourier  coefficients  w(k).  The  last  formula  shows  that  w(k )  are  determined  in 
terms  of  w(k),  by  folding  back  high  modes  on  the  lowest  ones,  due  to  the  discrete  resolution 
of  the  moments  of  w(x ):  all  modes  =  A[mod2/7  + 1]  are  aliased  at  the  same  place  since  they 
are  equal  on  the  gridpoints 


(1.2.12) 


D«(fc+p(2W+l))x^  _  :p(2W+l)r  _  «fcx„ 
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Let  us  rewrite  (1.2.7)  in  the  form 

w(k)  =  w{k )  +  ^  g«'p(2W+l)r  .  +  pf2 N  +  1)). 

Pr  0 

Returning  to  the  aliasing  error  in  (1.26),  we  now  have 


(1.2.13) 


Anw(x)  =  £ 
|fc|<AT 


53  e*p(2W+1)r  •  w(k  +  p  ■  (2N  +  1)) 


pr  0 


fcx 


We  note  that  Thw(x)  lies  outside  tcn  while  A^w(x)  lies  in  km,  hence  by  iP-orthogonality 
||ro(x)  —  ,ipNw(x)\\2H>  =  £  (1  +  l^!2)1  •  K*)|2  <—  truncation 

|fc|>  AT 

-f  £  (1  +  \k\2)a  ■  j  53  e,p{2*+1)r  ■  w(k  +  p(2N  +  1))|2  <—  aliasing. 

|Jfc|<w  o 

Both  contributions  involve  only  the  high  amplitudes  -  higher  than  N  in  absolute  value;  in 
fact  they  involve  precisely  all  of  these  high  amplitudes.  This  leads  us  to  aliasing  estimate 
£  (1  +  |A:|2)*|  £  e,p(2^+1)r  •  w(k  +  p{2N  +  1))|2  < 


|fc|<W 


p#o 


E  E(!  +  I*  +  p(2 N  +  1)|2)' l«(fc  +  p(2 If  +  l))!2- 

|fc|<JVp#0 


(1.2.15) 


■Maxw<Ar£ 

p*  o 


l  +  l*f 


\\Tnw(x)W2h,  •  £ 

p^  0 


l  +  |A:+p-(2W  +  l)|2 

r  i +n2  r 


< 


1  1-  4p2JV2 


Hence,  we  conclude  that  for  any  s  >  ~  we  have 
(1.2.16) 

Augmenting  this  with  our  previous  estimates  we  end  up  with  spectral  accuracy  as  before, 
namely 

(1.2.17) 

We  observe  that  ip^w(x)  is  nothing  but  the  trigonometric  interpolant  of  w(x)  at  the  equidis¬ 
tant  points  x  =  x 


MWx)!lff*  <  Const,,  •  ||7Vw(x)||h,,  s  >  ^ 


\w  -  ij}Nw\\H,  <  Const,  •  weH3[ 0,2tt),  s  >  a  >  i. 


N 


$NW{x)  |x=xm  =  £ 


(1.2.18) 


k=-N 


2N 


2 NTl^w(x-)e 


—ikx„ 


2W  i  N 

=  E  ■  577-7  E  =  «-(*,)• 

v=0  +  1  k=-N 
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This  shows  that  is  in  fact  a  t/>seudospectral  projection,  which  in  the  usual  sin-cos  formu¬ 
lation  reads 


(1.2.19) 


N 

ipNW  =  E)  'a-k  cos  kx  +  bk  sin  kx 
k= o 


2 

2N  +  1 


2  N 

E  wM 


i/=0 


cos  kxv 
sin  kxv 


Thus,  trigonometric  interpolation  provides  us  with  an  excellent  vehicle  to  perform  approxi¬ 
mate  discretizations  with  high  (=  spectral)  accuracy,  of  differential  and  integral  operations. 
These  can  be  easily  carried  out  in  the  Fourier  space  where  the  exponents  serve  as  eigenfunc¬ 
tion.  For  example,  suppose  we  are  given  the  equidistant  gridvalues,  wV)  of  an  underlying 
smooth  (i.e.,  also  periodic!)  function  w(x),w(x)eHa[Q,2n).  A  second-order  accurate  discrete 
derivative  is  provided  by  center  differencing 


Note  that  the  error  in  this  case  is,  0(h2)  =  w^(()h2 ,  no  matter  how  smooth  w(x)  is. 
Similarly,  fourth  order  approximation  is  given  via  Richardson  procedure  by 


dw(  \ 
S{x  =  I-)  = 


_  8[wu+1  -  uy.i]  -  [%;  -  w„-2’ 


12  h 


+  0(h4). 


The  pseudospectral  approximation  gives  us  an  alternative  procedure:  construct  the  trigono¬ 
metric  interpolant 


N 


(1.2.20) 


'ip^w(x)  =  E  w(k)e 
k=-N 


ikx 


2  N 


W 


w  = 


2  N  + 


T£ 

_ n 


wue 


—  ikxv 


u=0 


Differentiate  -  in  the  Fourier  space  this  amounts  to  simple  multiplication  since  the  exponen¬ 
tials  are  eigenfunctions  of  differentiation, 


d  N 

(1.2.21)  —'ipNw(x)  =  E  w(k)ike'kx, 

dx  k=-N 

and  we  approximate 

dw ,  d  . 

(1.2.22)  — (x  =  x„)  =  — t/>wu;(x)|I=Ii/  +  spectrally  small  error. 

dx  dx 

Indeed,  by  our  estimates  we  have  for  tu(x)ei-P[0,27r),s  >  1, 


(1.2.23) 


MaxcKiOwl^-u^z)  —  — V’wMaOl  —  Const.  |ju>(x)  —  iPnw(x)\\h^  < 


Const, 

N*-2 
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which  verifies  the  asserted  spectral  accuracy.  Similar  estimates  are  valid  for  higher  deriva¬ 
tives.  To  carry  out  the  above  recipe,  one  proceeds  as  follows:  starting  with  the  vector  of 
gridvalues,  w  =  (w0,  •  •  • ,  w2 jv),  one  computes  the  discrete  Fourier  coefficients 


i  2  N 

(1.2.24)  w(k)  =  2fi+1  I]  wve~ikx%  -N  <  k  <  N, 

or,  in  matrix  formulation 


(1.2.25) 


'  w(-N)  ' 

w0 

=  F 

: 

w(N) 

w2N  _ 

then  we  differentiate 


Fkv  = 


ikx„ . 


2N  +  1 


(1.2.26) 

or  in  matrix  formulation 
(1.2.27) 


w(k)  — >  ikw(k), 


'  w(-N)  ' 

'  w(-N)  ' 

'  -iN 

-»  A 

| 

,  A  = 

w(N) 

w(N ) 

iN 

and  finally,  we  return  to  the  “physical”  space,  calculating 


(1.2.28) 


or  in  matrix  formulation 


N 


E  ikw(k)e'kxu,  j/  =  0, 2 N, 


k=~N 


'  £<*•>  ' 

—iNw(—N) 

* 

=  F*  ■  (2N  +  1) 

• 

.  &*«> . 

iNw(N) 

(1.2.29) 

The  summary  of  these  three  steps  is 

(1.2.30) 


,  (2N  +  l)F*k  =  eikXv. 


w'(x0) 

w0 

l 

= 

: 

.  w'(x2N)  _ 

.  W2N  . 

t l>D  =  (2iV  +  1)F*AF, 


where  ipD  represents  the  discrete  differentiation  matrix,  and  similarly  i)Ds  for  higher  deriva¬ 
tives. 


Note:  Since  (2 N  -f-  1  )F*F  =  I2N+  i  (interpolation!)  we  apply  i>D3  =  (2 N  +  1)F*ASF.  How 
does  this  compare  with  finite  differences  and  finite-element  type  differencing? 
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In  periodic  second-order  differencing  we  have 


FD2  =  — 

2  2k 


0  1 
-1  0 

0 

1  0 


0 

-1 


fourth  order  differencing  yields 


FD* 


1 

m 


o 

-8 

1 


8  -1 
0 


-1 

8  -1 


-1 

0 

1 

0 


1  -8 
1 


-1 

0  8 

1  -8  0 


In  both  cases  the  second  and  fourth  order  differencing  takes  place  in  the  physical  space. 
The  corresponding  differencing  matrices  have  finite  bandwidth  and  this  reflects  the  fact  that 
these  differencing  methods  are  local.  Similarly,  finite-element  differencing, 


1  ,  4,1, 

r- + r- + 


Wu+i  -  Wv-1 

2  h 


corresponds  to  a  differencing  matrix 


FEi  = 


4 

6 

1 

6 

1 

6 


1 

6 


1 

6 

1 

6 

4 

6 


-l 


1 

2 h 


0  1 

-1  0 

0 

1  0 


0  -1 

0 


0 

-1 


1 

0 


We  still  operate  in  physical  space  with  O(N)  operations  (tridiagonal  solver)  and  locality 
is  reflected  by  a  very  rapid  (exponential  decay)  away  from  main  diagonal.  Nevertheless,  if 
we  increase  the  periodic  center  differences  stencil  to  its  limit  then  we  end  up  with  global 
pseudospectral  differentiation 

2N 

w^e 


^-1>Nw{xv)  =  Yj 

ax  k=-N 


(1.2.31) 

recall  the  Dirichlet  kernel  (1.1.33) 

(1.2.32) 


ik 


—  V  v,  *-'**<* 
.  2N  +  1  ~^Q 


Jkxv . 


N 

E  e' 

k=-N 


=  e 


e>(2N+i)x  _  sin(Ar  -f  |)x 


-  1 


sin  | 
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and  its  derivative, 


(l  2  as')  V  ike'kx  =  d  sm(N+  _  (N  +  I)  cos(N  +  |)z  sin  f  -  | cos  f  sin(N  +  \)x 
'  '  6  dx  sin  |  sin2 1 


(1.2.34) 


Np  ikeiH^)h  =  (N+l)cos[(N+\)(v-  ii)h] 

sin(^) 


Hence  (1.2.31),  (1.2.34)  give  us 


(1.2.35) 


„  ,  _  d  ,  ,  ,  ^  1  (-ir»  ,  ,  n,  (-1)*“" 

dx^wM  E2sin(f!qp.)  ^  !■"*  2sin(^)' 


In  this  case  if) D  is  a  full  (2iV  +  1)  x  (2W  +  1)  matrix  whose  multiplication  requires  0(N 2) 
operations;  however,  we  can  multiply  ^)D[u>]  efficiently  using  its  spectral  representation  from 
(1.2.30), 

iPD  =  (2iV  -f  l)F*AF. 

Multiplication  by  F  and  F*  can  be  carried  out  by  FFT  which  requires  only  0(N  log  N ) 
operating  and  hence  the  total  cost  here  is  almost  as  good  as  standard  “local”  methods,  and 
in  addition  we  maintain  spectral  accuracy. 

We  have  seen  how  the  pseudospectral  differentiation  works  in  the  physical  space.  Next, 
let’s  examine  how  the  standard  finite-difference/element  differencing  methods  operate  in 
the  Fourier  space.  Again,  the  essential  ingredient  is  that  exponentials  play  the  role  of 
eigenfunctions  for  this  type  of  differencing.  To  see  this,  consider  for  example  the  usual 
second  order  center  differencing,  D2(h),  for  which  we  have 


(1.2.36) 


D2{hy 


0ikxu+i  _  p—ikxv— i 


ism(kh)  ikx 


etKx\ 

c  \x=xl/) 


The  term  is  called  the  “symbol”  of  center  differencing.  By  superposition  we  obtain 

for  arbitrary  grid  function  (represented  here  by  its  trigonometric  interpolant) 


(1.2.37) 


i>Nw(x)  =  w(k)etk 


W./+1  ~  Wu-l 

2  h 


=  D2{h)ipNw  =  w(k)D2(h)e'kx\x=x 


(1.2.38) 


=  E 


r-  »sin(ik)  ,  , 


w{k)el 
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It  is  second-order  accurate  differencing  since  its  symbol  satisfies 


(1.2.39) 


ism(kh) 

h  =  * 


=  ik  +  0{k3h2). 


Note  that  for  the  low  modes  we  have  0(h2)  error  (the  less  significant  high  modes  are  differ¬ 
enced  with  0(1)  error  but  their  amplitudes  tend  rapidly  to  zero).  Thus  we  have 


(1.2.40) 


-  D2(h)i,Nwf  =  £ 

ax  k=-N  ^ 


<  Const. h*  ^  (1  +  jk|2)3|u'(k)|2  <  Const. hA  ■  H^ru/U^j, 

|ife|<W 


and  this  estimate  should  be  compared  with  the  usual 


%+l  - 


<  Const. A2  •  Max*  |u/3)(a:)| . 


The  main  difference  between  these  two  estiamtes  lies  in  the  fact  that  (1.2.41)  is  local,  i.e., 
we  need  the  smoothness  of  w{x )  only  in  the  neighborhood  of  x  =  xv  and  not  in  the  whole 
interval.  The  analogue  localization  in  the  Fourier  space  will  be  dealt  later. 

Similarly,  we  have  for  fourth  order  differencing  the  symbol 


.1  s'mkh  sin2&/t 

Z  —  4 - 

3  h  2  h 


=  ik  +  0{khhi). 


In  general,  we  have  difference  operators  whose  matrix  representation,  D,  is  of  the  form. 


(1.2.41) 


D  =  [dik]  —  N  <  j,k  <  N 


and  it  is  periodic  and  antisymmetric  (here  [£]  =  £[mod  2 N  -f  1]) 


(1.2.42) 


(i)  periodicity  djk  = 

(ii)  antisymmetry  djk  =  —  dkj. 


Matrices  satisfying  the  periodicity  property  are  called  circulant,  and  they  all  can  be  diago¬ 
nalized  by  the  unitary  Fourier  matrix 


(1.2.43) 


D  =  U*hU ,  U  =  (2iV  +  1)?  •  F,  U*U  =  I2N+1. 
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Indeed,  with  p  —  q  =  l  we  have 


N 


[U’DV]jk  =  -  E 


—ikxQ 


p,q=-N 


(1.2.44) 


ZiV  ^  1  As=-at 

2^  E  •  E 

"r  1  l,q=—N  q=-N 

(  0  j  T  k> 


=  < 


N 


E  «*%  j  =  K 

1  l=-N 


and  using  the  antisymmetry  we  end  up  with  symbols  A*. 


(1.2.45) 


N 


A  =  diag(A_w,  •  •  • ,  A# ),  Afc  =  2i  E  di  sin (klh). 

i=i 


As  an  example,  we  obtain  for  the  finite-element  differencing  system 


(1.2.46) 


x  . . /  4  1 

A‘  ~  ‘~(6  +  6e 


.sin  kh  /4  1 

(e  +  6 

6i  s'm(kh) 


tkh  +  -e~,kh  \  - 


h  4  +  2  cos  (kh) 


=  ik+  0(A4). 


In  general,  the  symbols  are  trigonometric  polynomials  or  rational  functions  in  the  “dual 
variable,”  kh,  which  has  “exact”  representation  on  the  grid  in  terms  of  translation  operator 
(polynomials  or  rational  functions),  and  accuracy  is  determined  by  the  ability  to  approximate 
the  exact  differentiation  symbol  ik  for|&|  ~  1,  see  Figure  2. 
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1.3.  Spectral  and  Pseudospectral  Approximations  -  Exponential  Accuracy 

We  have  seen  that  the  spectral  and  the  pseudospectral  approximations  enjoy  what  we 
called  “spectral  accuracy”  -  that  is,  the  convergence  rate  is  restricted  solely  by  the  global 
smoothness  of  the  data.  The  statement  about  “infinite”  order  of  accuracy  for  C°°  functions 
is  an  asymptotic  statement.  Here  we  show  that  in  the  analytic  case  the  error  decay  rate  is 
in  fact  exponential. 

To  this  end,  assume  that 

(1.3.1)  w(z)  =  ^2  w(k)e'kt,  | Im  z\  <  rj  <  r]0, 

Jess  —  oo 

is  27r-periodic  analytic  in  the  strip  —  t]q  <  Im  z  <  7j0.  The  error  decay  rate  in  both  the 
spectral  and  pseudospectral  cases  is  determined  by  the  decay  rate  of  the  Fourier  coefficients 
w(k).  Making  the  change  of  variables  (  =  elz  we  have  for 

(1.3.2)  v(C)  =  w(z  = 
the  power  series  expansion 

OO 

(1-3.3)  v(C)  =  £  w(k)(k. 

fc=-oo 


IT 


By  the  periodic  analyticity  of  w(z )  in  the  strip  \Imz\  <  rj  <  t]0,v(()  is  found  to  be  single¬ 
valued  analytic  in  the  corresponding  annulus 


(1.3.4) 


e  770  <  \tj\  <  e770, 


whose  Laurent  expansion  is  given  in  (1.3.3): 

«*>-5 


O~no 


<  r  <  e 


Vo 


This  yields  exponential  decay  of  the  Fourier  coefficients 

(1.3.6)  |w(*)|  <  M{rj)e-k\  M(V)  =  Max|Jmz|<>(.z)|,  0  <  V  <  Vo. 

We  note  that  the  inverse  implication  is  also  true;  namely  an  exponential  decay  like  (1.3.6) 
implies  the  analyticity  of  w(z).  Inserting  this  into  (1.1.24)  yields 


||u;  -  SatioII2  =  2tt  |ih(A:)|2  < 

|fc|>AT 


(1.3.7) 


<  2tt  •  M2(t))  •  y  e~2kr>  =  2tt^M  •  e~2N' 

z2r)  —  1 


\k\>N 


and  similarity  for  the  pseudospectral  approximation 


(1.3.8) 


w  —  -ip^w\\2  <  Const. 


M\n) 


0~2NV 


e2r)  —  1 

Note  that  in  either  case  the  exponential  factor  depends  on  the  distance  of  the  singularity 
(lack  of  analyticity)  from  the  real  line.  For  higher  derivatives  we  likewise  obtain 

,-2Nn 


(1.3.9) 


w  -  S^W |j^c  +  ||ty  -  TpNw\\2H<r  <  Const. N2<r  ■  M2(t]) 


e2n  —  1 


We  can  do  even  better,  taking  into  account  higher  derivatives,  e.g., 


(1.3.10) 

so  that  with 

(1.3.11) 

we  have 

(1.3.12) 
and  hence 

(1.3.13) 


kA{h)=hLJi^~ki^ 


i=o 


~kr) 


&|io(A:)j  <  Mi(Tj)e 


|w  -  5Vu>||//a  +  ||to  -  TpNwWn*  <  Const. M^tj) 


e~2Nr) 
e2r>  —  1 
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1.4.  The  Non-Periodic  Problem  -  Chebyshev  Approximation 
We  start  by  considering  the  second  order  SL  problem 


(1.4.1) 


—  Vl  —  x2—(Vl  —  —  1<X<1. 


This  is  a  special  case  of  the  general  SL  problem 

Noting  the  Green  identity 

(1.4.3)  (Lip,  (t>)u(x)  -  I  -(prpjtj)  +  =  p{x)[il> ,  <f>]\ba  +  {i>,  T^)w(x),  (V;>  =  i’fi  ~  W, 

J  a 

we  find  that  L  is  formally  self-adjoint  provided  certain  auxiliary  conditions  are  satisfied.  In 
the  nonsingular  case  where  p(a)  ■  p(b )  0,  we  augment  (1.4.2)  with  homogeneous  boundary 
conditions. 


(1.4.4) 


'ip(a)  =  <f>(a )  =  0,  ^(b)  =  =  0. 


Then  L  is  self-adjoint  in  this  case  with  a  complete  eigensystem  (A k,ipk(x)):  each 
w{x)eLu(x){a,  b]  has  the  “generalized51  Fourier  expansion 

(1.4.5)  w(x)  ~  p  w(k)ipk(x),  w{k)  = 


»(®)  ~  Y,  *(*0 


with  Fourier  coefficients 


(1.4.6) 


w(k )  = 


1  fb 


Wl*  Ja 


/  w(x)ipk(x)cj(x)dx. 
J  a 


The  decay  rate  of  the  coefficients  is  algebraic:  indeed 

™(k)  =  rrrr^  •  ~{L^k)w) u  = 

\\yk\\i  A k 


(1.4.7) 


=  WW  ’  +  ^k’  Lw ^ 


=  mi{p{x) '  £  ^  lU>w)& + ^  L(3)w)b 

The  asymptotic  behavior  of  the  eigenvalues  for  nonsingular  SL  problem  is 


Const  .k2 
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and  hence,  unless  w(x )  satisfies  infinite  set  of  boundary  restrictions  we  end  with  algebraic 
decay  of  w(k) 

1  p(x)  .  ,  ,l6  Const. 

This  leads  to  algebraic  convergence  of  the  corresponding  spectral  and  pseudospectral  pro¬ 
jections. 

In  contrast,  the  singular  case  is  characterized  by  having,  p(a)  =  p(b)  =  0,  in  which  case 
L  is  self-adjoint  independently  of  the  boundary  conditions  since  the  brackets  [  ,  ]  drop,  and 
we  have  spectral  decay-compare  (1.1.22) 

(L4'8) 

provided  w(x )  is  smooth  enough;  that  is,  the  decay  is  as  rapid  as  the  smoothness  of  w(x) 
may  permit. 

Returning  to  the  singular  SL  problem  (1.4.1)  we  use  the  transformation 


(1.4.9) 


did  Id 

x  =  cos  a.  —  =  -j—  •  —  - - =====—- 

dx  4§  dd  y/l  —  x2  d9 


which  yields 


(1.4.10) 


-  jjpW)  =  A  W)  =  ^(cos  6 ), 


obtaining  the  two  sets  of  eigensystems 


(1.4.11a) 


(A*  =  k2,  <f) k  =  cos  kd), 


(1.4.116) 


(A*.  =  k2,  (f>k  =  sin  kd). 


The  second  set  violates  the  boundedness  requirement  which  we  now  impose 


(1.4.12) 


|V,jfc(±l)l  <  Const., 


and  so  we  are  left  with 


(1.4.13) 


(A*  =  k2,Tpk(x)  =  cos^cos"1  x)). 


The  trigonometric  identity 


cos (k  +  l)d  =  2  cos  9  cos  kd  —  cos(&  —  1)0 
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yields  the  recurrence  relation 


(1.4.14) 


M+i(®)  =  2xi>k(x)  -  ipk-i{x),  i>o(x)  =  1 ,  V’i(a=)  =  z. 


hence,  ^.(x)  are  polynomials  of  degree  k  -  these  are  the  Chebyshev  polynomials 


(1.4.15) 


Tk(x)  =  cos(kcos~l  x ) 


which  are  orthonormal  w.r.t.  Chebyshev  weight  a>(x)  =  (1  —  x2)", 


(1.4.16) 


(ft«,  tax)).  =  £  ~0=r&  =  I  PX  =  \  i  =  k  >  o, 


j  7^  k, 


\\To\\l  =  7T  j  =  k  =  0. 


In  analogy  with  what  we  had  done  before,  we  consider  now  the  Chebyshev-Fourier  expansion 


(1.4.17) 


/  \  v-'  a  /  i\/n  /  \  a/l\  (M®)>  ^(x))^ 


To  get  rid  of  the  factor  |  for  &  =  0  we  may  also  write  this  as 


(1.4.18a) 


where 


(1.4.186) 


Ms)  ~  E  'M^)^(s) 

k=0 


(w(x),Tk(x))u  2  r1  w(x)cos(A:cos“1  x)dx 


Ikjx)),,  =  2  r* 

12  7 r  J-i 


'l  —  x2 


or  using  the  above  Chebyshev  transformation 


(1.4.18c) 


2  /■* 

!>(£)  =  —  w(cos  ^)  cos  &£  d(. 

7T  j£=0 


Thus,  we  go  from  the  interval  [—1,1]  into  the  27r-periodic  circle  by  even  extension,  with 
Fourier  expansion  of  w(cos  9),  compare  (1.1.9)  and  see  Figure  3, 

1  f2*  2  f* 

w(k)  =  —  ty(cos^)cos  k(d£  =  —  w(cos  £)  cos  k(d(. 

7T  0  7 r  J(=  0 

Another  way  of  writing  this  employs  a  symmetric  doubly  infinite  Fourier-like  summation, 
where 


(1.4.19) 


1  00 

Ms)  ~  n  E  ™(k)Tk{X) 
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w(9 )  =  w(x  =  cos9 ) 


Figure  3: 


with  T_ jt(x)  =  T/t(x)  and 
(1.4.20)  w( 


2  f1  w(x)Tk(x) 

(k)  —  —  >  -  K-^-dx,  —  oo  <  k  <  cc. 

W  W-i 


The  Parseval  identity  reflects  the  completeness  of  this  system 


l*(o)l2  +  ^EK*)l’ 

=  JEK*)I2 

*  fc=0 

which  yields  the  error  estimate 

II™  -  S/HIt  =  tE  l™(*OI2- 

4  fc>AT 

Now,  in  order  to  get  a  measure  of  the  spectral  convergence  we  have  to  estimate  the  decay 
rate  of  Chebyshev  coefficients  in  terms  of  the  smoothness  of  w(x)  and  its  derivatives;  to 
this  end  we  need  Sobolev  like  norms.  Unlike  the  Fourier  case,  {Tfc(x)}  is  not  complete  with 
respect  to  Ha  -  orthogonality  is  lost  because  of  the  Chebyshev  weight.  So  we  can  proceed 
formally  as  before,  see  (1.1.24), 

(1.4.22)  \\w  -  =  2x  £  |*(fc)|2  <  £ 

k>N  k>N  \L  +  iV  ) 

i.e.,  if  we  define  the  Chebyshev-Sobolev  norm 

IMIbj.  =  E(1  +  l*l2)*W*)l2, 

fc=0 


11  t  MI2  fw2(x)dx  1 

ll™(*)Hr  =  /  7f==f  =  4 


(1.4.21) 
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then  we  have  spectral  accuracy 


||u>  —  >Sjv"U>||:r  <  Const,  •  — ,  weH^ [— 1, 1]. 

N3 

In  fact  the  H?  space  can  be  derived  from  appropriate  inner  product  in  the  real  space  as  done 
in  Fourier  expansion.  The  correct  inner  product  is  given  by  -  compare  (1.1.19) 

3  3  r2*  d2p  d2p 

(1.4.23)  (wl,w2)H2.  =  Y^(Lpm,Lpw2)T  ^  Y  je=0J§^w^cos9)j^w^cos6)d6' 

p=°  x=co  S0P=°  a(' 


so  that 


(1.4.24) 


0  j  k, 

(^fc>  Tj)H2T>  =  <  tj.  • 

—  Y  ^4p>  j  —  k  (with  7 r  factor  at  j  =  k  =  0). 
.  2p= o 


Hence  the  Fourier  coefficients  in  this  Hilbert  space  behave  like 


(1.4.25) 


and  the  corresponding  norm  is  equivalent  to 


(1.4.26) 


Mlk*  ~  Y(l  +  fc2)25K*)ls 


The  reason  for  the  squared  factors  here  is  due  to  the  fact  that  L  is  a  second  order  differential 
operator,  unlike  the  first-order  D  =  ^  in  the  Fourier  case,  i.e., 


(1.4.27) 


£(i  +  ~  £ 

k= 0  p=0 


involves  the  first  2s-derivatives  of  w{x)  -  appropriately  weighted  by  Chebyshev  weight.  This 
completes  the  analogy  with  the  Fourier  case,  and  enables  us  to  estimate  derivative  as  well- 
compare  (1.1.28  -  1.1.29), 


(1.4.28) 


|u;  -  <  Const,-— a  <  s,  weH3T [-1, 1]. 


Next,  let’s  discuss  the  discrete  setup.  Because  we  seek  an  even  extension  of  the  upper  semi¬ 
circle  we  consider  the  case  of  even  number  of  grid  points  -  equally  distributed  along  the  unit 
circle.  One  choice  is  to  look  at 


_  ,  /  J. »  71  /  ,  n  u, 

6v  =  r  +  vh  =  {v+-)-  h=-,r  =  - 
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which  gives  the  usual  Chebyshev  interpolation  points,  see  Figure  4.  The  second  choice  - 
which  we  will  concentrate  from  now  on  -  takes  into  account  the  boundaries,  where  we  look 
at 

,  7T  .  2tT  . 

e„  =  «h  =  (k  =  — ,  r  =  0) 

which  yields,  see  Figure  5, 

(1.4.29)  xv  —  cos(i/h),  v  =  0, 1,  •  •  •  ,N. 


The  trapezoidal  rule  is  replaced  by  Gauss-Lobatto  rule  for  (1.4.20),  i.e.,  starting  with 
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(1.4.18c),  we  have 


(1.4.30)  w(k )  =  —  f  w(cos£)  cos  k£d£  »  —  V)  "u^cos  6V)  cos  , 

Tt  J(=  0  TT  IV 

and  we  end  up  with  the  discrete  Chebyshev  coefficients 


(1.4.31) 


2  ^ 

=  0  <  k  <  N. 

^  i/=0 


Indeed,  this  corresponds  to  the  case  of  Fourier  interpolation  with  even  number  of  equidistant 
gridpoints  9V  in  (A.  1.2),  for 

i  iN  ,  n  r  ,  0_ 

•«*-“**  =  +  ^ 

Z7r  i/=0  71  v=0 

2  w 

=  —  E  "w»  cos(k6„). 

^  i/=0 

Then  one  may  construct  the  Chebyshev  interpolant  at  these  N  +  1  gridpoints 


(1.4.32) 


ipNw(x)  =  J2  "w(k)Tk(x). 


We  have  an  identical  aliasing  relation 


(1.4.33) 


w(k)  =  ^  w{k  +  2pN) 


(compare  A. 1.5),  and  hence  spectral  convergence,  i.e.,  (compare  (1.2.14)  and  (1.2.12)) 


(1.4.34) 


|u;(x)  —  iPnw(x)\\h°  <  Const,  •  — — ,  weH^y  s  >  a, 


where  Const,  ~  ||^||h' • 


Example:  We  have  the  Sobolev  imbedding  of  L°°  space 


m*)i < o  e  \m\  <  dDi+wrE 


(1  4-  k 2)° 


<  Const,,  •  ||ttf||s?>  a  >  \. 


Consequently, 


Maxx|u;(x)  -  ipNw(x) |  <  (Const,  ~  |M|h*)  •  3  >  a  > 
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In  particular,  with  s  =  N  +  1  we  obtain  the  ’’usual”  estimate  for  the  near  minmax  approxi¬ 
mation  (  collocated  at  at  xu  =  cos  77)) 

e~N 

Maxx|zu(a;)  —  f>Nw(x)\  <  Const. ||ty||H^+i  •  — y. 

We  briefly  mention  the  exponential  convergence  in  the  analytic  case.  To  this  end  we  employ 
Bernstein’s  regularity  ellipse,  Er,  with  foci  ±1  and  sum  of  its  semi  axis  =  r.  Denoting 


(1.4.35) 


M(rj)  =  Max^mz)!,  r  =  en 


We  have 


Theorem.  Assume  w(x)  is  analytic  in  [-1,1  J  with  regularity  ellipse  whose  sum  of  semiaxis 
=  r0  =  e™  >  1.  Then 

\\w(x)  —  ^iu(aj)||#<,  -{- ||-u>(2;)  —  Snw(x) \\2h<,  <  Const.  ^ |  ■  N2cre~2N’1. 


Proof:  The  transformation  2  =  ((  +  (  1)/2  takes  Ero  from  the  z-plane  into  the  annulus 
ro1  <  ICI  <  ro  in  the  C-plane.  Hence,  u(()  =  2w  {z  =  admits  the  power  expansion 

(1.4.36)  v(()  =  2w  =  £  w(k)(k,  V 1  <  |d  <  r0  =  eCo; 

indeed,  setting  (  =  el8  and  recalling  w(—k )  =  w(k),  the  above  expansion  clearly  describes 
the  real  interval  [-1,1] 

OO 

(1.4.37)  w(z  —  cos  fl)  =  ■ w(k )  cos  k6. 

k=0 

By  Laurent  expansion  in  (1.4.36) 

(1.4.38)  m  =  -2 
hence 

(1.4.39)  |ui(A)|  <  M( v)e~kTt 

and  the  result  follows  along  the  lines  of  (1.3. 7-8). 

Finally,  we  conclude  with  discussion  on  Chebyshev  differencing.  Starting  with  grid  values 
wv  at  Chebyshev  points  xv  =  cos  (^77))  one  constructs  the  Chebyshev  interpolant 

N  9  N 

(1.4.40)  ipNw(x)  =  "w(k )  ■  Tk(x),  w(k)  =  —  ^  "wu  cos (kxu). 

k=o  v=0 
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One  can  compute  w(k),  0  <  k  <  N,  efficiently  via  the  cos-FFT  with  0(N  log  N)  operations. 
Next,  we  differentiate  in  Chebyshev  space 


(1.4.41) 

In  this  case,  however,  Tk{x)  is  not  an  eigenfunction  of  instead  ~Tk{x)  -  being  a  polyno¬ 
mial  of  degree  <  k  —  1,  can  be  expressed  as  linear  combination  of  {Tj(a;)}^Zo  0n  fact  2*(x) 
is  even/odd  for  even/odd  k's):  with  c0  =  2,  cj.>0  =  1  we  obtain 


(1.4.42) 

and  hence 

(1.4.43) 

Rearranging  we  get 

d 


r/^  =  7  E  W), 

dx  Cfc  o<j<k 

k—j  odd 


^Niu{x)  =  Y™{k)—  Y  tAx)- 

k=  0  Cfc  o<]<k 

k—j  odd 


dx 


N 


(1.4.44) 


~1pNw(X)  =  Y  "™'{k)Tk(x),  #'(*)  =  —  Y  p*(p) 

ax  k- o  Ck  p>k+ i 

p+fc  odd 


and  similarly  for  the  second  derivative 


(1.4.45) 


*"(*)  =  —  Y  p(p2  ~  k2)Hp)- 

ck  p>k+2 
p+k  even 


The  amount  of  work  to  carry  out  the  differentiation  in  this  form  is  0(N2)  operations  which 
destroys  the  NlogN  efficiency.  Instead,  we  can  employ  the  recursion  relation  which  follows 
directly  from  (1.4.44) 

(1.4.46)  w'(k  +  1)  =  w\k  —  1)  •  Ck-i  —  2kw(k). 

To  see  this  in  a  different  way  we  note  that 

sin(&  -f  1)#  =  sin (k  —  1)9  -f  2  sin  9  cos  k8, 


which  leads  to 


1  dTk+i  1 


k  +  1  dx  k  —  1  dx 


+  2Tk(x), 


and  hence 
d 


dx 


ip?rw(x)  = 


f:"kw(k)ln(x) = 

k=o  K 

N 


1  ^  1 

-  Y  "{'u)'{k  —  1)  —  w'(k  +  l))-Tl.(x )  =<■ —  summation  by  parts 
2  Jb=0  k 

N  N 


=  i  £  "2 ib'(k)Th(x)  =  £  "w\k)Tt(x) 


k= 0 


k=0 
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as  asserted.  In  general  we  have 


(1.4.47) 


w^(k  4-1)  =  w^{k  -  !)<:*_!  -  2kw^-1)(k). 


With  this  w{k)  can  be  evaluated  using  0(N)  operations,  and  the  differentiated  polynomial 
at  the  grid  points  is  computed  using  another  cos-FFT  employing  0(N  log  N)  operations 


(1.4.48) 


d  N 

—TpNw(x)\x=Xl/  =  ]T"tI/(fc)  cos  kxv 

ax  k=o 


with  spectral/exponential  error 
(1.4.49) 


d  d 

MaxI=IJ—  w(x)  -  —  i})Nw(x)\  <  ^ 


Const,  •  ~  ^  <  o  <  s, 

Const,  •  e~^r] 


The  matrix  representation  of  Chebyshev  differentiation,  DT,  takes  the  almost  antisymmetric 
form 

f  <=n-i),+‘ 


(DT  ),k  = 


Ck  Xj  -  xk 
X 


3  l  k, 


~W^T)  ;  =  M(0,iV)’ 


2N2  4- 1 

6 

21V2  4-1 


j  =  k  =  0, 
j  =  k  =  N. 
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2.  TIME  DEPENDENT  PROBLEMS 


2.1.  Initial  Value  Problems  of  Hyperbolic  Type 
The  wave  equation, 

(2.1.1)  wtt  =  a2wXX) 

is  the  prototype  of  P.D.E.’s  of  hyperbolic  type.  We  study  the  pure  initial-value  problem 
associated  with  (2.1.1),  augmented  with  27r-periodic  boundary  conditions  and  subject  to 
prescribed  initial  conditions, 

(2.1.2)  w(x,0)-f(x),  u)t(x,0)  =  g(x). 


We  can  solve  this  equation  using  the  method  of  characteristics,  which  yields 

f (x  -f  at)  +  / (x  —  at)  1  rx+at 


(2.1.3) 


w 


(x,t)  = 


+ 


2a 


rx^at 

/  g(s)ds. 

J  x— at 


We  shall  study  the  manner  in  which  the  solution  depends  on  the  initial  data.  In  this  context 
the  following  features  are  of  importance. 


1.  Linearity:  the  principle  of  superposition  holds. 

2.  Finite  speed  of  propagation:  influence  propagates  with  speed  <  a.  This  is  the  essential 
feature  of  hyperbolicity.  In  the  wave  equation  it  is  reflected  by  the  fact  that  the  value  of 
w  at  ( x ,  t)  is  not  influenced  by  initial  values  outside  domain  of  dependence  ( x—at ,  x+at ) 

3.  Existence  for  large  enough  set  of  addmissible  initial  data:  arbitrary  (7^°  initial  data 
can  be  prescribed  and  the  corresponding  solution  is  (7£°. 

4.  Uniqueness:  the  solution  is  uniquely  determined  for  — oo  <  t  <  co  by  its  initial  data. 

5.  Conservation  of  Energy.  The  wave  equation  (2.1.1)  describes  the  motion  of  a  string 
with  kinetic  energy,  ~p  fw2dx,  and  potential  one,  |T  /  w2xdx,  (T / p  =  a2).  In  order 
to  show  that  the  total  energy 

^Totai  =\pj  (wt  + 

is  conserved  in  time  we  may  proceed  in  one  of  two  ways:  either  by  the  so  called  energy 
method  or  by  the  Fourier  method. 
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The  Energy  Method. 


Rewrite  (2.1.1)  as  a  first  order  system 


(2.1.4a) 


d_ 

dt 


Ui 

'  0 

a 2 

d 

Ui 

Ui 

y-2 

1 

0  _ 

dx 

u2  _ 

> 

U2 

dw 

dt 

dw 

dx 


or  equivalently, 

(2.1.46) 

The  essential  ingredient  here  is  the  existence  of  a  positive  symmetrizer,  H  >  0, 


du  du 
dt  dx 


HA  = 


0  a2 
a2  0 


=  A,  —  Aj ,  H  = 


1  0 
0  a 2 


(2.1.5) 

so  that  multiplication  by  H  on  the  left  gives 

(2.1.6)  Hut  =  A,ux. 

Now,  multiplying  by  ur  we  are  led  to 

(2.1.7)  (u,Hut)  =  (u,Atux), 

and  the  real  part  of  both  sides  are  in  fact  perfect  derivatives,  for  by  the  symmetry  of  H, 

Re(u,  Hut )  =  i(u,  Hut )  +  ~{Hut,  u )  = 


tt  \  1,"  rr  N  9 

-(u,Hut)+  ~(ut,  Hu)  =  — 
and  similarily,  by  the  symmetry  of/1,,  we  have 


1  1  d  ‘  1 

Re(u,  A,ux)  =  -(u,/!,^)^-  ~{A,ux,u )  =  —  ~(u,A,u) 

Hence,  by  integration  over  the  27r-period  we  end  up  with  energy  conservation,  asserting 


(2.1.8)  —  /  (wf  +  a2wl)dx  =  —  f  (u,Hu)dx  =  /  ~(u,  A,u)dx  =  0. 

dt  Jx  dt  Jx  Jx  Ox 


d 


We  note  that  the  positi  vi  iy  of  wets  not  used  in  the  proof  and  is  assummed  just  for  the 
sake  of  making  (u,  Hu)  an  admissible  convex  “energy  norm.” 
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The  Fourier  Method. 


Fourier  transform  (2.1.4b)  to  get  the  O.D.E. 
(2.1.9) 

whose  solution  is 


fill 

— (fc,i)  =  ikAu(k,t), 


(2.1.10) 


uik,  t ) 


-  JkAtz. 


u(k ,  0), 


where  u(k,  0)  is  the  Fourier  transform  of  the  initial  data.  Now,  for 


(2.1.11) 

we  find 
(2.1.12) 
or 

(2.1.13) 


A  =  TAT 


-l 


A  = 


—a 

,T  = 

—a  a 

a 

j 

7 

1  1 

u(M)  =  Te^T-^M) 


T~lu(k,t ) 


,-t  kat 
0 


T-'tyk,  0) 


and  hence  (since  the  diagonal  matrix  on  the  right  is  clearly  unitary),  the  £2-norm  of 
T~lu(k,t )  is  conserved  in  time,i.e., 


(2.1.14)  ||r-1i(M)||!  =  ||r-1i(t,o)||2,  r-1  =  -1 


1  —a 
-1  —a 


Summing  over  all  modes  and  using  Parseval  we  end  up  with  energy  conservation 

[,  2  ,  2  2U  4  2  [  {Wt~awx\2  fwt  +  awx\ 2 

+  0  “’■)*  =  4“  J,  +  (-rjH  dx 

=  4o2  J  ||r-'n||2<ix  =  87m!  Y)  j|7— 'e(A,  i)f  =  Const. 


(2.1.15a) 


as  asserted. 

We  note  that  the  only  tool  used  in  the  energy  method  was  the  existence  of  a  positive  sym- 
metrizer  for  A,  while  the  only  tool  used  in  the  Fourier  method  was  the  real  diagonalization 
of  A]  in  fact  the  two  are  related,  for  if  A  =  TAT~l  then  for  H  =  {T~x)*T~l  >  0  we  have 


(2.1.156) 


HA  =  (T  1yAT~1  =  As  =  A^)  Areal  diagonal. 


Energy  conservation  implies  (in  view  of  linearity)  uniqueness,  and  serves  as  a  basic  tool 
to  prove  existence.  It  will  be  taken  as  definition  of  hyperbolicity.  It  implies  and  is  implied 
by  the  qualitative  properties  (1)  -  (4). 
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We  now  turn  to  consider  general  P.D.E.’s  of  the  form 

A  d  Q 

(2.1.16a)  =  P(x,t,D)u,  P(x,t,D)  =  Y,Aj(x,t)—, 

o— i  J 

with  27r-periodic  boundary  conditions  and  subject  to  prescribed  initial  conditions, 

(2.1.166)  u(x,0)  =  f(x). 

We  say  that  the  system  (2.1.16a)  is  hypberbolic  if  the  following  a  priori  energy  estimate 
holds: 

(2.1.17)  IKM)I|l2(x)  <  Const. j  •  |(u(a;,  0)||^2(a:),  -T  <t<T. 

As  we  shall  see  later  on  this  is  equivalent  to  energy  conservation  with  appropriate  energy 
renorming. 

Here  are  the  basic  facts  concerning  such  systems. 

The  Constant  Coefficients  Case. 


(2.1.18) 


=  P(D)u,  P(D)  =  ^2Aj—,  Aj  =  constant  matrices. 
j= l  c)x 


dt 


Define  the  Fourier  symbol  associated  with  P(D ): 

d 

(2.1.19)  P(ik)  =  i  J2  Ajkjt  k  =  (kt,  k2,  •  •  • ,  kd)eRd, 

j=i 

which  naturally  arises  as  we  Fourier  transform  (2.1.18), 

ps 

(2.1.20)  —u(k,t)  =  P(ik)u(k,i). 

Sovling  the  O.D.E.  (2.1.20)  we  find,  as  before,  that  hyperbolicity  amounts  to 

(2.2.21)  ]|e^(‘fc)t||  <  Constj,  -T<t<T,  for  ally's. 

For  this  to  be  true  the  necessary  Garding-Petrovski  condition  should  hold,  namely 


(2.1.22) 


JReA[P(zA:)]|  <  Const. 


Example:  For  the  wave  equation,  A[.P(i£)]  =  ±ika. 

But  the  Garding-Petrovski  condition  is  not  sufficient  for  the  hyperbolic  estimate  (2.1.17)  as 
told  by  the  counterexample 


Ui 

a 

1  ' 

d 

Ui 

. u 2  . 

0 

a 

dx 

.  U2  . 
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As  before,  in  this  case  we  have  A[P(iA:)]  =  ±ika,  hence  the  Garding-Petrovski  condition  is 
fulfilled.  Yet,  Fourier  analysis  shows  that  we  need  both  |[?x1(x,  0)Jj^2(x)  and  |j4j^(a:,0)||i,2(3C)  in 
order  to  upperbound  jju^x,  t)  j|L2(*)-  Thus,  the  best  we  can  hope  for  with  this  counterexample 
is  an  a  priori  estmate  of  the  form 

IK®,  OIUj(x)  <  ConstT  •  ||u(®,  0)||Hi(l),  -T  <t<T. 

We  note  that  in  this  case  we  have  a  ’’loss”  of  one  derivative.  This  brings  us  to  the  notion  of 
weak  hyperbolicity. 

We  say  that  the  system  (2.1.16a)  is  weakly  hyperbolic  if  there  exists  an  s  >  0  such  that  the 
following  a  priori  estimate  holds: 

IK®,  t)  ||ij(l)  <  Constr  •  ||«(x,  0)[|h-(i),  -T  <t<T. 

The  Garding-Petrovski  condition  is  necessary  and  sufficient  for  the  system  (2.1.18)  to  be 
weakly  hyperbolic.  The  necessary  and  sufficient  characterization  of  hyperbolic  systems  is 
provided  by  the  Kreiss  matrix  theorem.  It  states  that  (2.1.21)  holds  iff  there  exists  a  positive 
symmetrizer  H(k)  such  that 

(2.1.23)  R e[H(k)P(ik)}  =  0,  0  <  m  <  H(k)  <  M, 

and  this  yields  the  conservation  for  j|rz(x,  ^)il//  =  2»El!!*(*)ili,„.  i-e.. 

2n  Y^{u(k,  t)),  H(k)u(k,  t )) 

k 

is  conserved  in  time. 

Remark:  For  an  a  priori  estimate  forward  in  time  (0  <  t  <  T ),  it  will  suffice  to  have 

(2.1.24)  R e[H(k)P(ik)}  =  ^[H(k)P(ik)  +  P(ik)H(k))  <  0. 

Indeed,  we  have  in  this  case 

li(u(k),  H(k)u(k))  <  (Re[7J (k)p(ik)}u(k),u(k))  <  0. 
and  hence  summing  over  all  k’s  and  using  Parseval’s 

Special  important  cases  are  the  strictly  hyperbolic  systems  where  P(ik)  has  distinct  real 
eigenvalues,  so  that  P(ik)  can  be  real  diagonalized 

P(k )  =  fT(fc)A(ifc)T-1(A:), 
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and,  as  before,  H(k)  =  ( T~1(k))*T~1(k )  will  do.  The  other  important  case  consists  of 
symmetric  hyperbolic  systems  which  can  be  symmetrized  in  the  physical  space,  i.e.  there 
exists  an  H  >  0  such  that 

HAj  =  Aj.  =  A%. 

Most  of  the  physically  relevant  systems  fall  into  these  categories. 


Example:  Shallow  water  equations  (linearized) 


u 

,  d 

u 

,  d 
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Vo  . 

can  be  symmetrized  with 


H  = 


The  Variable  Coefficient  Case. 


(2.1.25)  -^  =  P(x,t,D)u. 

This  is  the  motivation  for  introducing  the  notion  of  hyperbolicity  as  is:  freeze  the  coefficients 
and  assume  hyperbolicity  of  ut  —  P(x0)t0,  D)u  uniformly  for  each  (x0,t0) ;  then  (unlike  the 
case  of  weak  hyperbolicity),  the  variable  coefficients  problem  is  also  hyperbolic. 

Remark:  This  result  is  based  in  the  invariance  of  the  notion  of  hyperbolicity  under  low- 
under  perturbations;  it  restricts  hyperbolic  system  to  be  of  first-order. 

So  far  we  have  dealt  with  hyperbolicity  via  the  Fourier  method  studying  the  algebraic 
properties  of  its  symbol  P(ik)\  we  can  also  work  with  the  energy  method. 


For  example,  if  we  assume  that  P(x,t,  D)  is  semi-bounded,  i.e.,  if 
(2.1.26)  -M|MiL(l)  <  Re(«,  P(x,  t,  D)u)L2{x)  <  M||u||^(l),  0  <  M, 

then  we  have  hyperbolicity. 
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Example:  The  symmetric  hyperbolic  case  Aj(x,t)  =  Aj(x,t):  we  can  rewrite  such  sym¬ 
metric  problems  in  the  equivalent  form 


du 

¥ 


^  £;( AiU) 


+  Bu, 


dAj 
dxj ' 


In  this  case  the  symmetry  of  the  Aj's  implies  that  |  +  Sy  is  skew- 

adjoint,  i.e.,  integration  by  parts  gives 


=  0. 

L2(x) 


Therefore  we  have 

Re(ti,  jP(x,  t,  Z))ii)^,2(x)  =  Re(-8ti,  u)£2(x)> 

and  hence  the  semi-boundedness  requirement  (2.1.26)  holds  with  M  =  ||ReJ5|| .  Conse¬ 
quently,  if  Aj(x,t )  are  symmetric  (or  at  least  symmetrizable)  then  the  system  (2.1.16a)  is 
hyperbolic. 
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2.2.  Initial  Value  Problems  of  Parabolic  Tipe. 


The  heat  equation, 

(2.2.1)  ut  =  auxx,  a  >  0, 

is  the  prototype  of  P.D.E.’s  of  parabolic  type.  We  study  the  pure  initial-value  problem 
associated  with  (2.2.1),  augmented  with  27T-periodic  boundary  conditions  and  subject  to 
initial  conditions 

(2.2.2)  u(x,0)  =  f(x). 

We  can  solve  this  equation  using  the  Fourier  method  which  gives 

(2.2.3)  u(M)  =  e~ak7tf(k). 

It  shows  the  dissipation  effect  (=  the  rapid  decay  of  the  amplitudes,  |w(A:,  i)|,  as  functions  of 
the  high  wavenumbers,  j /cj  1)  in  this  case,  which  is  the  essential  feature  of  parabolicity. 

As  before,  we  study  the  manner  in  which  the  solution  depends  on  its  initial  data. 

1.  Linearity:  the  principal  of  superposition  holds. 

2.  Uniqueness:  the  solution  is  uniquely  determined  for  t  >  0  by  the  explicit  formula 

1  -i2 

(2.2.4)  u(x,t)=  Q(x  —  y,t)f(y)dy,  Q{z)  =  -j=e  ««*  >  0. 

J yzz—oc  V  4tt at 

3.  Existence  for  large  enough  set  of  addmisible  initial  data:  bounded  initial  data  f(x) 
can  be  prescribed  (or  at  least  j/(x)|  <  eux*),  and  the  corresponding  solution  is  C°°  - 
in  fact  u(x,  t  >  0)  is  analytic  because  of  exponential  decay  in  Fourier  space. 

4.  The  maximum  principle:  follows  directly  from  the  representation  of  u(x ,  t )  as  a  convo¬ 
lution  of  f(x)  with  the  unit  mass  positive  kernel  Q(z). 

5.  Energy  decay:  as  usual  we  may  proceed  in  one  of  two  ways. 

The  Fourier  Method.  We  start  with 

ll|)«(x,  Oil,  <2*E  l/MI’ '  ' '  |e-*’f  ]  <  Const.r*  •  ||/||2fa, 

(since  the  quantity  inside  the  above  brackets  is  maximized  at  k2at  =  s ).  The  last  a  priori 
estimate  shows  that  the  parabolic  solution  becomes  infinetly  smoother  than  its  initial  data 
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is  as  we  ’’gain”  infintely  many  s-derivetives,  and  at  the  same  time,  higher  derivetives  decay 
faster  as  t  f  oo.  Alternatively,  we  can  work  with  the 

Energy  Method. 


Multiply  (2.2.1)  by  u  and  integrate  to  get 


(2-2.6)  25l|tt,kM  - 

and  in  general 

/n  n  „\  1  d  ,.dsu..0  ,.d3+1u..0 

(2-2-7)  j  sfc111*  -  'Const'Ha?TTll^ 

successive  integration  of  (2.2.7)  yields  (2.2.5). 

Turning  to  general  case,  we  consdier  mth-order  P.D.E.’s  of  the  form, 


1^11^112  /  n _ <.  ii^4+1uii2 


(2.2.8) 


=  P(x,t,D)u,  P(x,t,D)  =  yy  A3{X)t)D3. 


We  say  that  the  system  (2.2.8)  is  weakly  parabolic  of  order  a  if 


(2.2.9) 


-u(X|f)l|i2  <  Const.r|i|/ajju(a;,0)||L2(l). 


In  the  constant  coefficients  case  this  leads  to  the  Garding-Petrovski  characterization  of 
parabolicity  of  order  (3 ,  requiring 


ReA  P(ik)  =  y  Aj(ik)J  <  -A  ■  \kf  +  B. 
bl=o 

Remark:  Generically  we  have  a  ~  (3  =  m  the  order  of  dissipation  which  is  necessarily 


The  extension  to  the  variable  coefficients  case  (with  Lipschitz  continuous  coeeficients) 
may  proceed  in  one  of  two  ways.  Either,  we  freeze  the  coeeficients  and  apply  the  Fourier 
method  to  the  constant  coeficients  problems;  or,  we  may  use  the  energy  method,  e.g.,  inte¬ 
gration  by  parts  shows  that  for 

F('X’ =  +B,ey  +  Cu’ 

with  A]  +  A)  >  8  >  0,  and  Bj  =  B*,  the  corresponding  systems  (2.2.8)  is  parabolic  of 
order  2. 

Example:  ut  =  auxx  +  uxxx  is  weakly  parabolic  of  order  two,  yet  does  not  satisfy  Petrovski 
parabolicity. 
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2.3.  Well-Posed  Time-Dependent  Problems. 


Hyperbolic  and  parabolic  equations  are  the  two  most  important  examples  of  time-dependent 
problems  whose  evolution  process  is  well-posed.  Thus,  consider  the  initial  value  problem 


(2.3.1) 


du  .  . 

—  =  P(x,t,D)u. 


We  assume  that  a  large  enough  class  of  admissible  initial  data 


(2.3.2) 


u(x,t  =  0)  =  f(x) 


there  exists  a  unique  solution,  u (cc ,  i).  This  defines  a  solution  operator,  E(t,r)  which  de¬ 
scribes  the  evolution  of  the  problem 


(2.3.3) 


u(t)  =  E(t,  r)u(r). 


Hoping  to  compute  such  solutions,  we  need  that  the  solutions  will  depend  continuously  in 
their  initial  data,  i.e., 


(2.3.4) 


Ju(t)  -  u(t)jj  <  ConstrjJu(O)  -  u(0)||h*  0  <  t  <  T. 


In  view  of  linearity,  this  amounts  to  having  the  a  priori  estimate  (boundedness) 


(2.3.5) 


\u(t)  =  E(t,  t)u(t) ||  <  Constr|ju(r)|]H*,  0  <  t  <  T, 


which  includes  the  hyperbolic  and  parabolic  cases. 

Counterexample:  (Hadamard)  By  Cauchy-Kowaleski,  the  system 


du  du 

-  +  A-  =  0,  u  = 


,A  = 


0  +1 

-1  0 


has  a  unique  solution  for  arbitrary  analytic  data,  at  least  for  sufficiently  small  time.  Yet, 
with  initial  data 


(2.3.6) 


we  obtain  the  solution 


,  s  inni 

ui(x,0)  = - ,  u2(x,0)-0, 

n 


,  ,,  cosh  nt  sin  nx  sinhntcosnx 

ux(x,i)= - ,  u2{x,t)  = - 

n  n 

which  tends  to  infinity  ||w(-,<)||n_00  — »  oo,  while  the  initial  data  tend  to  zero.  Thus  this 
Laplace  equation,  —  0,  is  not  well-posed  as  an  initial- value  problem. 

Finally,  we  note  that  a  well-posed  problem  is  stable  against  perturbations  of  inhomogeneous 


38 


data  in  view  of  the  following 

Duhamel’s  principle.  The  solution  of  the  inhomogeneous  problem 


(2.3.8)  —  =  P(x,  t,  D)u  +  F(x,  t) 
is  given  by 

(2.3.9)  u(t)  =  E(t,  OHO)  +  J*  E(t,  r)F(r)dT. 

Proof:  We  have 

=  £w.0M0)]  +  |[£o£(<.r)F(iMT] 

=  P(x,i,D)\E(t,0)u(0)]  +  E(t,t)F(t)  +  j‘_^[E(t,T)F(t)]dr 

=  P(x,t,  D)[E{t,  0)u(0)  +  f^  E(t,T)F(T)ir\  -v  F(t)  =  P(x,t,D)u(t)  +  F. 

This  implies  the  a  priori  estimate 

(2.3.10)  IHOII  <  Constr||u(0)||  +  Constr  /  \\F(t)\\Hs(It,  0  <  t  <  T, 

Jr= 0 

as  asserted. 
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3.  THE  FOURIER  METHOD  FOR  HYPERBOLIC  AND  PARABOLIC  EQUA¬ 
TIONS 


3.1.  The  Spectral  Approximation 


We  begin  with  the  simplest  hyperbolic  equation  -  the  scalar  constant-coefficients  wave 
equation 


(3.1.1) 


du  du 
dt  a  dx 


subject  to  initial  conditions 


(3.1.2) 


u(x,0)  =  f(x), 


and  periodic  boundary  conditions. 

This  Cauchy  problem  can  be  solved  by  the  Fourier  method:  with  f(x )  =  f(k)etkx 
we  obtain  after  integration  of  (3.1.1), 


(3,1.3) 

with  solution 


( k ,  t)  =  ikau(k,  t), 


(3.1.4)  ii(k,t)  =  eikatf(k), 

and  hence 


(3.1.5)  u(x,t)  —  ^2  e'hatf(k)eikx  =  £  f(k)eik^  =  /(*  +  at). 

k  k 

Thus  the  solution  operation  in  this  case  amounts  to  a  simple  translation 


E(t,  t)u( x,  r)  =  u(x  +  a(t  -  r),  t),  || E(t,  r)||  =  1. 


This  is  reflected  in  the  Fourier  space,  see  (3.1.4),  where  each  of  the  Fourier  coefficients  has 
the  same  change  in  phase  and  no  change  in  amplitude;  in  particular,  therefore,  we  have  the 
a  priori  energy  bound  (conservation) 


(3.1.6)  |W'.f)l|2  =  2*  £  I Hk,  t) I2  =  2,  £  \m2  =  ll/(-)H2. 

k  k 


We  want  to  solve  this  equation  by  the  spectral  Fourier  method.  To  this  end  we  shall  ap¬ 
proximate  the  spectral  Fourier  projection  of  the  exact  solution  Sun  =  Snu(x,  t).  Projecting 
the  equation  (3.1.1)  into  the  N- space  we  have 


(3.1.7) 


8un 

dt 


=  S, 


N 


du 

% 
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Since  5 at  commute  with  multiplication  by  a  constant  and  with  differentiation  we  can  write 
this  as 


(3.1.8) 


dupj  dux 
dt  a  dx 


Thus  Un  =  Snu  satisfies  the  same  equation  the  exact  solution  does,  subject  to  the  approxi¬ 
mate  initial  data 


(3.1.9)  u,v(*  =  0)  =  Sat/. 

The  resulting  equations  amount  to  2 N  -f  1  ordinary  differential  equations  (O.D.E.)  for  the 
amplitudes  of  the  projected  solution 

(3.1.10)  —u^{k,  t)  =  ikau^(k,  t ),  —  N  <  k  <  N, 

subject  to  initial  conditions 

(3.1.11)  uN(k,0)  =  f(k). 

Since  these  equations  are  independent  of  each  other,  we  can  solve  them  directly,  obtaining 

(3.1.12)  uN(k>t)  =  e,katf(k) 
and  our  approximate  solution  takes  the  form 

(3.1.13)  M*,t)=  £)  f(k)e^x+at\ 

k=-N 

Hence,  the  approximate  solution  u^(x)  t)  =  /a?(x  -f  at)  satisfies 

(3.1.14)  u( x,  t)  -  uk(x,  t )  =  E(t,  0 )f(x)  -  E(t ,  O)Stff(x) 


and  therefore,  it  converges  spectrally  to  the  exact  solution,  compares  (1.1.26), 


(3.1.15) 


IK*)  ~  «jv(0ll  ^  IW>0)(7  -  5'jv)/(x)||  < 

<  IK  -  SN)f(x)\\  <  Const||/||H.  •  ~ • 


Similar  estimates  holds  for  higher  Sobolev  norms;  in  fact  if  the  initial  data  is  analytic  then 
the  convergence  rate  is  exponential.  In  this  case  the  only  source  of  error  comes  from  the 
initial  data,  that  is  we  have  the  error  equation 


(3.1.16) 


d  d 

— [u  -  ittf]  =  a— [u  -  Utf] 


v\ 
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subject  to  initial  error 


(3.1.17) 


U  —  Un{t  =  0 )  —  f  —  Jn- 


Consequently,  the  a  priori  estimate  of  this  constant  coefficient  wave  equation 


u  -  <  ConstT|j/  -  fN ||  <  Const. ||/||/f3  ■  —  Constr  =  1. 


(3.1.18) 

Now  let  us  turn  to  the  scalar  equation  with  variable  coefficients 

(3.1.19) 

This  hyperbolic  equation  is  well-posed:  by  the  energy  method  we  have 


du  du  . 

—  =  a(x,  <)— ,  a(x,  t)  =  Z7r  —  periodic. 


=  - /  uza(x,t)u-J  a-(x,t)u'- 

l  f  A  -V  ^ 

(3.1.20)  J^u2(x,t)dx  =  J  ua(x,t)ux  =  -- J  ax(x,  t)u2(x,  t)dx, 

and  hence 


(3.1.22a)  IK®,  0IU2(x)  <  Constr  •  ||/(s)|| 

with 

(3.1.226)  Constj  =  eMT,  M  =  Maxx>t[— ax(x,  <)]. 

In  other  words,  we  have  for  the  solution  operator 


||5(«,T)u(r)||i2(x)  <  eM^1  r)||u(r)||Ml) 

and  similarly  for  higher  norms.  As  before,  we  want  to  solve  this  equation  by  the  spectral 
Fourier  method.  We  consider  the  spectral  Fourier  projection  of  the  exact  solution  uN  = 
S^u(x,t)\  projecting  the  equation  (3.1.19)  we  get 

(3.1.23)  frtUN  ~ 

Unlike  the  previous  constant  coefficients  case,  now  SN  does  not  commute  with  multipli¬ 
cation  by  a(x,  t),  that  is,  for  arbitrary  smooth  function  p(x,  t)  we  have  (suppressing  time 
dependence) 

N  /  oo 

(3.1.24a)  SNa(x)p(x)  =  a(fc-j)p(j) 

k=:—N  \j=—oo 


.  .du 
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while 


£  a (k-j)p(j) 

k— —co  \j=—N 

Thus,  if  we  exchange  the  order  of  operations  we  arrive  at 
(3.1.25)  =  a(x,t)~~  -  [a(x,t)SN  -  SNa(x,t )]^. 

While  the  second  term  on  the  right  is  not  zero,  this  commutater  between  multiplication  and 
Fourier  projection  is  spectrally  small,  i.e., 

||5'ivo(x)/t)(a;)  -  a(x)5^(*)||£J(«)  = 


(3.1.245)  cl(x)Snp(x )  =  £ 


(3.1.26) 


||(SV  -  I)a(x)p(x)  +  a(x)(7  -  5w>(x)||i2(x)  < 

1 

<  Const. ||a(x)p(x)||/f.  •  jjj  +  Const. ||a(a;)||ioo(l)  •  ||p(x)|| 


1 


and  so  we  intend  to  neglect  this  spectrally  small  contribution  and  to  set  as  an  approximate 
model  equation  for  the  Fourier  projection  of  u(x,  t) 


Ovn 

dt 


a(x,  t) 


OVN 

dx 


Yet,  the  second  term  may  lie  outside  the  N-space,  and  so  we  need  to  project  it  back,  thus 
arriving  at  our  final  form  for  the  spectral  Fourier  approximation  of  (3.1.19) 


(3.1.27) 


dvN 

dt 


Again,  we  commit  here  a  spectrally  small  deviation  from  the  previous  model,  for 

(3.1.28)  1|(7  -  Stf)ap(x)j|L2(x)  <  Const  j|a(x)p(x)|j,/.  •  — . 

The  Fourier  projection  of  the  exact  solution  does  not  satisfy  (3.1.22),  but  rather  a  near  by 
equation, 

(3.1.29)  ^^  =  SN[a(x,t)^~^^rFN{x,t) 


where  the  truncation  error.  Fn(x,  t )  is  given  by 


(3.1.30) 


Fu{x,t)  —  Stf 


a{x,  t)(I  -  Sn) 


du 

dx 


43 


The  truncation  error  is  the  amount  by  which  the  (projection  of)  the  exact  solution  misses 
our  approximate  mode  (3.1.27);  in  this  case  it  is  spectrally  small  by  the  errors  committed  in 
(3.1.26)  and  (3.1.18).  More  precisely  we  have 

(3.1.31)  \\Fn(x,  i)||Lj(x)  ]|a(x,  t)||i,2(x)  •  IM|h»+1 

depending  on  the  degree  of  smoothness  of  the  exact  solution.  We  note  that  by  hyperbolicity, 
the  later  is  exactly  the  degree  of  smoothness  of  the  initial  data,  i.e.,  by  the  hyperbolic 
differential  energy  estimate 

(3.1.32)  li-^Xz,  *)IU»(*)  <  IK*,  0IUs(*) '  ll/lk*+1  •  Jj~ 

and  in  the  particular  case  of  analytic  initial  data,  the  truncation  error  is  exponentially  small. 

From  this  point  of  view,  our  spectral  approximation  (3.1.27)  satisfies  an  evolution  model 
which  is  spectrally  away  from  that  of  the  Fourier  projection  of  the  exact  solution  (3.1.29). 
This  is  in  addition  to  the  spectrally  small  error  we  commit  initially,  as  we  had  before 

(3.1.33)  Vflr(t  =  0)  =  Sxf  =  /jv- 

We  now  raise  the  question  of  convergence.  That  is,  whether  the  accumulation  of  spectrally 
small  errors  while  integrating  (3.1.27)  rather  than  (3-1.29),  give  rise  to  an  approximate 
solution  v^(xy  t)  which  is  only  spectrally  away  from  the  exact  projection  u^(x,  t).  We  already 
know  that  the  distance  between  u^(x,  t)  and  the  exact  solution  u(x,  t)  -  due  to  the  spectrally 
small  initial  error  -  is  spectrally  small  as  we  have  seen  in  the  previous  constant  coefficient 
case. 

To  answer  this  convergence  question  we  have  to  require  the  stability  of  the  approximate 
model  (3.1.27).  That  is,  we  say  that  the  approximation  (3.1.27)  is  stable  if  it  satisfies  an  a 
priori  energy  estimate  analogous  to  the  one  we  have  for  the  differential  equation 

(3.1.34)  ||ujv(0ll  <  Const.eMij|vAr(0)||. 


Clearly,  such  a  stability  estimate  is  necessary  in  any  computational  model.  Otherwise,  the 
evolution  model  does  not  depend  continuously  on  the  (initial)  data,  and  small  rounding 
errors  can  render  the  computed  solution  useless.  And  on  the  positive  side  we  will  show  that 
the  stability  implies  the  spectral  convergence  of  an  approximate  solution  ujv(x,f).3  Indeed 
the  error  equation  for  e^(t)  =  u^(t)  -  v^(t)  takes  the  form 


(3.1.35) 


dtN  c 

1T  =  Sn 


den 

dx 


+  FN(x,t) 


3We  note  that  in  the  previous  constant  coefficient  case,  the  approximate  model  coincides  with  the  differ¬ 
ential  case,  hence  the  stability  estimate  was  nothing  but  the  a  priori  estimate  for  the  differential  equation 
itself. 
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Let  E^{t,  r )  denote  the  evolution  operator  solution  associated  with  our  approximate  model. 
By  the  stability  estimate  (3.1.34) 

(3.1.36)  ||£n(*,tMt)||  <  Conste"^>|Mr)||. 

Hence,  by  (3.1.36)  together  with  Duhammel’s  principle  we  get  for  the  inhomogeneous  error 
equation  (3.1.35) 

(3.1.37a)  eN(t)  =  EN(t,  0)eN(0)  +  [  EN(t,T)FN{r)dT 

J-r=0 

and 

(3.1.376)  ||ejv(*)||  <  Const.eMt  ||e^(0)|| l2(x)+  I  || FN(xyT)\\L^x)dT  . 

Jt=0 

In  our  case  e^(0)  =  /at  —  S/m  =  0,  and  the  truncation  error  F^(x,t)  is  spectrally  small; 
hence 

(3.1.38)  ||e//  =  itAr(f)  —  fAr(t)||  <  Const. eMt  ■  — 

where  the  constant  depends  on  j|a(m, i)t|z,oo(0  anc*  |!/I!hj+1>  i-e-i  restricted  solely  by  the 
smoothness  of  the  data.  In  the  particular  case  of  analytic  data  we  have  exponential  conver¬ 
gence 

(3.1.39)  |M<)  =  uN(t)  -  v*(i)||  <  Const. eWt  •  e^N. 


Adding  to  this  the  error  between  u^(t)  and  u(i )  (-  which  is  due  to  the  spectrally  small  error 
in  the  initial  data  between  f#  and  /)  we  end  up  with 

f  ^7  for  Hs+l  initial  data 

(3.1.40)  |ju(t)  —  'yN'(^)l!  ^  Const. eMt  •  < 

(  e-r,w  for  analytic  initial  data 


To  summarize,  we  have  shown  that  our  spectral  Fourier  approximation  converges  spectrally 
to  the  exact  solution,  provided  the  approximation  (3.1.27)  is  stable. 

Is  the  approximation  (3.1.27)  stable?  That  is,  do  we  have  the  a  priori  estimate  (3.1.34)? 
To  show  this  we  try  to  follow  the  steps  that  lead  to  the  analogue  estimate  in  the  differential 
case,  compare  (3.1.20).  Thus,  we  multiply  (3.1.27)  by  vn{x,t)  and  integrate  over  the  in¬ 
period,  obtaining 


(3.1.41) 
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But  Vtf(x,t)  is  orthogonal  to  (/  —  Sn)  so  adding  this  to  the  right-hand  side  of 

(3.1.41)  we  arrive  at 

(3.1.42)  =  J^vN(x,t)a(x,t)^-dx 

and  we  continue  precisely  as  before  to  conclude,  similarly  to  (3.1.22),  that  the  stability 
estimate  (3.1.34)  holds 

(3.1.43)  ||vjv(£)il  <  Const.eMt||uw(0)||,  M  =  MaxX)t[-oa.(s, i)]. 


In  the  constant  coefficient  case  the  Fourier  method  amount  to  a  system  of  (2N  +  1) 
uncoupled  O.D.E.’s  for  the  Fourier  coefficients  of  v ^  =  u jy  which  were  integrated  explicitly. 
Let’s  see  what  is  the  case  with  problems  having  variable  coefficients  say,  for  simplicity, 
a  =  a(x).  Fourier  transform  (3.1.22)  we  obtain  for  v(k,t)  =  v^^kyt)  -  the  kth-Fourier 
coefficient  of  v^(x,t)  =  Y^k=-N  t)e'kx, 


(3.1.44) 


-- =  Ea(^  iYiKh  t)>  -N  <k<N. 
at 


In  this  case  we  have  a  (2N  +  1)  x  (2 N  +  1)  coupled  system  of  O.D.E.’s  written  in  the 
matrix- vector  form,  consult  (1.2.46) 

j  \  v{~N,t)  1 


(3.1.45)  — =  AAv(t)}  v(t)  = 

at 


Akj  =  d(k  -  j),  A  =  diag(zk). 


We  can  solve  this  system  explicitly  (since  a  (•)  was  assumed  not  to  depend  on  time) 


(3.1.46) 


i >(t)  =  eAAtf>(0); 


that  is,  we  obtain  an  explicit  representation  of  the  solution  operator 


(3.1.47) 


EN{t,T)  =  Fj  eA^FN,  A  =  An,A.  =  A; 


where  F A  denote  the  spectral  Fourier  projection 

'  v(-N )  1 

(3.1.48)  Fnvn(x)  =  • 

v(N) 

We  note  that  in  view  of  Parseval  identity  ||i'Vtw(x)||2  =  ||u/^(s)||£2(x)  (modulo  factorization 
factor),  hence,  stability  amounts  to  having  the  a  priori  estimate  on  the  discrete  symbol 
EN(t,r)  =  eANh^~T\  requiring 

(3.1.49)  ||e^A(t-T)j|  <  Const.e^"^. 
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The  essential  point  of  stability  here,  lies  in  having  a  uniform  bound  in  the  RHS  of  (3.1.49) 
independent  on  the  order  of  the  system,  e.g.,  a  straightforward  estimeate  of  the  form 


(3.1.50) 


||e^ArA(f—T) ||  < 


will  not  do  because  j^Zco00'  The  essence  °f  the  a  priori  estimate  we  obtained  in 
(3.1.22),  and  likewise  in  (3.1.42),  was  that  the  (unbounded)  operator  P(x,t,D)  =  a(x,t)dx 
is  semi-bounded,  i.e., 

Re  a(s,i)^-  =~  -  ~-(a(xtty)  =  ~ax(x,t)] 


namely,  compare  (2.1.26) 


^Re  a(x,  £)—  <  M\\u\\l2(l 

VI  OXj  J  L2(x) 


and  likewise  for  R e  (5V  ja(:r,  f)  In  the  present  form  this  is  expressed  by  the  sharper 
estimate  of  the  matrix  exponent,4  compare  (3.1.50) 


(3.1.51) 


jeiwA(t-r)jj  <  gUReAjvAlKt-r) 


This  time,  |!ReAjvA|]  like  the  Re[P(a;,  t,  D)],  is  bounded.  Indeed,  [ReAA]^  =  | \a(k  —  j)ij  -f 
a(j  —  k)ik],  and  since  a(x,t )  is  real  (hyperbolicity!)  then  a(p)  =  a(—p),  i.e., 


(3.1.52) 


[ReAA]^  =  -i(j  -  k)a(k  -  j)  -N<j,k<  N. 


Thus,  ReAA  is  a  Toeplitz  matrix,  namely  its  (k,j)  entry  depends  solely  on  its  distance  from 
the  main  diagonal  k  —  j;  we  leave  it  as  an  exercise  -  using  our  previous  study  on  circulent 
matrices  in  (1.2.43)  -  to  see  that  its  norm  does  not  exceed  the  sum  of  absolute  values  along 
the,  say,  zeroth  (j  =  0)  row,  i.e., 


(3.1.53) 


IIReiwAII  <  r  E  M*)l 

1  k=  AT 


which  is  bounded,  uniformly  with  respect  to  N,  provided  a(x,  t)  is  sufficiently  smooth,  e.g., 
we  can  take  the  exponent  M  to  be 


(3.1.54) 


m  =  -  E  M*)l  <  ji/E^WW-  E  tt< 

L  k=-N  z  v  k=-N  K 


<  ||a*,(®,  0IUa(x) 


4To  see  this,  use  Duhammel’s  principle  for  =  ReAA5(i)  +  F(l)  where  F(t)  =  iImAAeAM  or  integrate 
directly. 
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which  is  only  slightly  worse  than  what  we  obtained  in  (3.1.43). 

A  similar  analysis  shows  the  convergence  of  the  spectral-Fourier  method  for  hyperbolic 
systems.  For  example,  consider  the  N  x  N  symmetric  hyperbolic  problem 


(3.1.55a) 


d  du 

—  =  A(x,  t)~  +  B( x,  t)u ,  with  symmetric  A( x,  t). 

C/L  C/JO 


We  note  that  if  the  system  is  not  in  this  symmetric  form,  then  (in  the  1-D  case)  we  can  bring 
it  to  the  symmetric  form  by  a  change  of  variables,  i.e.,  the  existence  of  a  smooth  symmetric 
H(x,  t )  such  that  H(x ,  t)A( x,  t )  is  symmetric,  implies  that  for  w(x,  t)  —  T_1(x,  t)u(x ,  t)  with 
H  =  (2’-1)*J’-1  we  have,  compare  (2.1.15b) 


(3.1.555) 


^  =  T~1(x,t)A(x,t)T(x,t)^  -f  C(x,t)w(x,t) 


where  T~1(x,t)A(x,t)T(x,t)  =  T*(x,t)H(x,t)A(x,  t)T(x,t)  is  symmetric,  and  C(x,t )  = 
B(x,  t)+^J-(x,t)  —  T~1(x,  t)A(x,  t)^(x,t).  The  spectral  Fourier  approximation  of  (3.1.55a) 
takes  the  form 


=  SN  A( x,  t)~z—  +  SNB(x,  t)vN(x,  t). 


(3.1.56) 


Its  stability  follows  from  integration  by  parts,  for  by  orthogonality 

(3.1.57a)  J  v2N(x,t)dx  =  J  v^A(x,t)^^-dx  +  J  u^B(x,t)u^dx  <  M  J  v2N(x,t)da 


where 


(3.1.575) 


and  hence 


(3.1.58) 


M  =  Max,,,  -  +  ReB(*.‘) 


IMOIU’W  <  ="‘IMo)||. 


The  approximation  (3.1.56)  is  spectrally  accurate  with  (3.1.55)  and  hence  spectral  conver¬ 
gence  follows.  The  solution  of  (3.1.56)  is  carried  out  in  the  Fourier  space,  and  takes  the 
form 


(3.1.59) 


d  N  - 

-7-v(k,t)  =  Y  A(k  -N  <k<  N, 

ai  j=-N 


which  form  a  coupled  (2 N  +  1)  x  (2 N  +  1)  system  of  O.D.E.’s  for  the  ( 21V  +  l)-vectors  of 
Fourier  coefficients  v(k,  t). 

There  are  two  difficulties  in  carrying  out  the  calculation  with  the  spectral  Fourier  method. 
First,  is  the  time  integration  of  (3.1.59);  even  in  the  constant  coefficient  case,  it  requires  to 
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compute  the  exponent  eA'At  which  is  expensive,  and  in  the  time-dependent  case  we  must 
appeal  to  approximate  numerical  methods  for  time  integration.  Second,  to  compute  the 
RHS  of  (3.1.59)  we  need  to  multiply  an  (2 N  +  1)  x  (2 N  +  1)  matrix,  AA  by  the  Fourier 
coefficient  vector  which  requires  0{N'i)  operations.  Indeed,  since  A  is  a  Toeplitz  matrix  and 
A  is  diagonal,  we  can  still  carry  out  this  multiplication  efficiently,  i.e.,  using  two  FFT’s  which 
requires  0(N  log  N)  operations.  Yet,  it  still  necessitates  to  carry  out  the  calculation  in  the 
Fourier  space.  We  can  overcome  the  last  difficulty  with  the  pseudospectral  Fourier  method. 

Before  leaving  the  spectral  method,  we  note  that  its  spectral  convergence  equally  applies 
to  any  P.D.E. 

du 

(3.1.60)  =  P(x,t,D)u 

with  semi-bounded  operator  P(x,  t,  D ),  e.g.,  the  symmetric  hyperbolic  as  well  as  the  parabolic 
operators.  Indeed,  the  spectral  approximation  of  (3.1.60)  reads 

(3.1.61)  J^.  =  sNP{x,t,D)vN. 

Multiply  by  vN  and  integrate  -  by  orthogonality  and  semi-boundedness  we  have 

(3-1.62)  \jtJx vlf(xit)dx  =  ^(vN.P(x,t,D)vN)  <  M  J^v2N(x,t)dx. 

Hence  stability  follows  and  the  method  converges  spectrally. 
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3.2.  The  Pseudospectral  Approximation 


We  return  to  the  scalar  constant  coefficient  case 


(3.2.1) 


du  du 
dt  a  dx 


subject  to  periodic  boundary  conditions  and  prescribed  initial  data 


(3.2.2) 


«(*»0)  =  f(x). 


To  solve  this  problem  by  the  pseudospectral  Fourier  method,  we  proceed  as  before,  this  time 
projecting  (3.2.1)  with  the  pseudospectral  projection  to  obtain  for  u ^  —  ip^u(x,t) 

,  .  dun  ,  (  du\ 

(3-2.3)  -3T  =  hu  • 


Here,  ipfj  commutes  with  multiplication  by  a  constant,  but  unlike  the  spectral  case,  it  does 
not  commute  with  differentiation,  i.e.,  by  the  aliasing  relation  (1.2.2)  we  have 

(3.2.4a)  E(M*))e'“  =  E  E*!*  +  +  l))p{k+j(2N  +  l))eifa 


where  as 


(3.2.46) 


k=-N  j 


E  ikY,P\k  +  H2N  +  \)]eib. 

OX  k=-N  k=-N  j 


The  difference  between  these  two  operators  is  a  pure  aliasing  error,  i.e.,  we  have  for  = 
Spj  +  Apr,  see  (1.2.13) 

-  j-(^Arp)  =  An,4~  P  =  Y  Yilk  +  K2N +  l))p(k  +  j{2N +  l)]elkx 
U'<"  C‘  ‘  L  X  J  k=-N  j^tO 

which  is  spectrally  small.  Sacrificing  such  spectrally  small  errors,  we  are  led  to  the  pseu¬ 
dospectral  approximation  of  (3.2.1) 


(3.2.5) 

subject  to  initial  conditions 

(3.2.6) 


8vn  dvpj 

~W  =  a~dT 


Vfi, r(t  =  0)  =  ijjMf- 


Here,  v x  =  u//(a :,t)  is  an  N-degree  trigonometric  polynomial  which  satisfies  a  nearby  equa¬ 
tion  satisfied  by  the  interpolant  of  the  exact  solution  ipNu(x,t).  That  is,  u at  =  ^n(a;,f) 
satisfies  (3.2.5)  modulo  spectrally  small  truncation  error 


(3.2.7) 


~W  ~  aYx~  +  FN^X,t^  FN{x^)~a^N 
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where  by  (3.2.3),  F#(x,t)  =  a\ipN  g  -  £(i}>Nuj\,  and  by  (1.2.17)  it  is  indeed  spectrally 
small 


(3.2.8) 


||-FMM)II  <  M  ~  ^n)u}\\}  <  |aj  ||u||H.+i~ 


The  stability  proof  of  (3.2.5)  follows  along  the  lines  of  the  spectral  stability,  and  spectral 
convergence  follows  using  Duhammel’s  principle  for  the  stable  numerical  solution  operator. 
That  is,  the  error  equation  for  e #  =  u n  —  is 


deu  de^  ,  „  ,  ,N 

-ar=“-& +F"(:M) 


(3.2.9) 
whose  solution  is 

(3.2.10)  eN(t)  -  EN(t,0)(fN  -  i>Nf)  +  [  EN{t,r)FN{x,r)dT. 

Jt= 0 

Hence,  by  stability 

(3.2.11)  |Mi)||  <  Const.  eMt  •  \\u\\H.«  ~  <  Const.eMt||/||H,+l  • 
this  together  with  the  estimate  of  the  pseudospectral  projection  yields 


(3.2.12) 


jju(t)  —  vN(t)\\  <  Const. e 


Mt 


_L 

N‘ 


for  7P+1  initial  data 


e~r)N  fQr  analytic  initial  data 


To  carry  out  the  calculation  of  (3.2.5)  we  can  compute  the  discrete  Fourier  coefficients  v(k,  t) 
which  obey  the  O.D.E., 


(3.2.13) 


dv 

dt 


(k,v)  =  ikav(k ,  t), 


as  was  done  with  the  spectral  case;  alternatively,  we  can  realize  our  approximate  interpolant 
vw(x,t)  at  the  2 N  4- 1  equidistant  points  x„  —  uh,  and  (3.2.5)  amounts  to  a  coupled  (2 N+  1) 
-  O.D.E.  system  in  the  real  space 


(3.2.14) 


-^(Xv,  t)  =  a~^(x  =  xv,  t)  v  =  0,1,  —  ,2iV. 


(3.2.15) 


xw(x»,o)  —  y('^i')- 


Let  us  turn  to  the  variable  coefficient  case, 

du 


(3.2.16) 


.  du 

dt  =  a(x'l)Tx 
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The  pseudospectral  approximation  takes  the  form 


(3.2.17a) 


dvN  .  ,,dvN 

-gr  =  ^w  a(X,t)-%r 


(3.2.17  6) 


vN(xu,0)  =  /(s„). 


It  can  be  solved  as  a  coupled  O.D.E.  system  in  the  Fourier  space,  but  just  as  easily  can  be 
realized  at  the  2 N  +  1  so-called  collocation  points 


(3.2.18a) 


dvN(xv,  t)  _  t)<k#(x  =  x  t) 

dt  -  L)  l) 


(3.2.186) 


vN{xv,t  =  0)  =  /(®„). 


The  truncation  error  of  this  model  is  spectrally  small  in  the  sense  that  Ujv  =  ipifU  satisfies 


(3.2.19) 


where 


(3.2.20) 


9un  ,  \  <  n 

~q^  =  Fn  a(x,t)-^-  +  FN(xit) 


du"  d 

Fn(x,  t)  =  i>N  a(x,t)—  -ipN  a(x,t)—(ipNu) 


is  spectrally  small 


(3.2.21) 


l|i?w(x,0il  <  \\i>N  a(*,<)^-(/‘”^)u]l!  - 


<  ||a(x>t)||L»  •  ||/|ftf. 


Hence,  if  the  approximation  (3.2.12)  stable,  spectral  convergence  follows.  Is  the  approxima¬ 
tion  (3.2.12)  stable?  The  presence  of  aliasing  errors  is  responsible  to  a  considerable  difficulty 
in  proving  this,  and  currently  this  is  an  open  question.  If  we  try  to  follow  the  differential 
and  spectral  recipe,  we  should  multiply  by  vu(x,t)  and  integrate  by  parts.  However,  here 
vx(x,t)  is  not  orthogonal  to  (/  -  V’iv)[- '  •]  which  otherwise  would  enable  us  to  follow  the 
differential  estimate  of  f  v^(x,  t)a(x,  t)^-(x,  t)dx  in  terms  of  fxvj:(x,t)dx;  more  precisely, 
we  have  for  ipu  =  +  An  that  I  —  ipu  =  I  —  Sn  —  An  where  /  vn(I  —  5tv)[-  •  • }dx  =  0, 

compare  (3.1.41),  (3.1.42);  yet  J  vnAn[-  •  -]dx  leaves  us  with  an  additional  contribution  which 
cannot  be  bounded  in  terms  of  fxvlf(x,t)dx  in  order  to  end  up  with  the  stability  proof  with 
Gronwall’s  inequality.  To  shed  a  different  light  on  this  difficulty,  wre  can  turn  to  the  Fourier 
space;  we  write  (3.2.17)  in  the  form 

,  .  dv M  ,  . dvM 


(3.2.22) 


dvif  dvif 

ar  =  a(x)-sT 
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and  Fourier  transform  to  get  for  the  kth  Fourier  coefficient 


(3.2.23a) 


N 

E  a{k  —  j,  t)ijv(j,  t) 

j=-N 


i.e., 

(3.2.236)  £v(t)  =  ANAv{t)  Akj  =  £  a[k  -  j  +  p(2N  +  1)]. 

at  „ 


This  time,  Re^A  is  unbounded.  This  difficulty  appears  when  we  confine  ourselves  to  the 
discrete  framework:  multiplying  (3.2.18a)  by  v(xv,t)  and  trying  to  sum  by  parts  we  arrive 
at 


(3.2.24) 


=  E^  l-a{x,t)v\x,t) 


X—Xv 


y  2a  (a'^>  t)vN(Xv>  Os 


but  the  first  term  on  the  right  does  not  vanish  in  this  case  -  it  equals,  by  the  aliasing  relation, 
to 


E 


d  ri 


dx 


-a(x,t)v  (x,t) 


-I 


^H  +  £p(2tf  +  l)^[p-(2tf  +  l)] 

xr=xv  0  . 


and  a  loss  of  one  derivative  is  reflected  by  the  factor  2W+1  inside  the  right  summation.  There 
are  two  main  approaches  to  enforce  stability  at  this  point:  skew-symmetric  differencing  and 
smoothing.  We  discuss  these  approaches  in  the  next  two  subsections. 
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3.3  Skew-Symmetric  Differencing  The  essential  argument  of  well-posedness  for  symmetric 


hyperbolic  s)rstems  with  constant  coefficients  is  the  fact  that  (say  in  the  1-D  case)  P(D)  = 

A-^  is  a  skew-adjoint  operator  what  is  loosely  called  “integration  by  parts” .  With  variable 

coefficients  this  is  also  true,  modulo  low-order  bounded  terms,  i.e. , 

q  1  3  3  "  1 

(3.3.1)  P(x,t,D)  =  A(x,t)—  =  -  A(x,t)—A—{a(x,t)-)  -  -Ax(x,t). 

The  stability  proofs  of  spectral  methods  follow  the  same  line,  i.e.,  we  have  in  the  Fourier 
space,  compare  (3.1.45),  (3.1.52) 

(3.3.2)  An  A  =  -  [AatA  —  +  -  |v4atA  +  AMatJ 

and  stability  amounts  to  show  that  the  second  term  in  (3.3.2)  is  bounded:  for  then  we  have 
in  (3.3.2)  (as  in  (3.3.1))  a  skew-adjoint  term  with  an  additional  bounded  operator.  The 
difficulty  with  the  stability  of  pseudo-spectral  methods  arises  from  the  fact  that  the  second 
term  on  the  right  of  (3.3.2)  is  unbounded, 

(3.3.3)  limAr_ool!^(-divA  -f  A*i;V)|!  T  oo. 

To  overcome  this  difficulty,  we  can  discretize  the  symmetric  hyperbolic  system  (again,  say 
the  1-D  case) 

,  .  du  A.  .du 

(3.5.4a)  -  =  A(x>t)- 

when  the  spatial  operator  is  already  put  in  the  “right”  skew-adjoint  form,  compare  (3.3.1), 
(3.3.46)  %==\  A(X,t^  +  ”  ^Ax(x,t)u. 

The  pseudospectral  approximation  takes  the  form 

3  I  (  3  3  "l 

(3.3.5)  •  =  -  i’w  +  Q^^N{A(x,t)uN)  -  ~^N{Ax{x)t)uN). 

In  the  Fourier  space,  this  gives  us 

dv  1  1  8An  . 

(3.3.6)  Tl  - 

Now,  An  A  +  A  An  is  symmetric  because  A  is,  js  bounded  and  stability  follows. 
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3.4.  Smoothing 


We  have  already  met  the  process  of  smoothing  in  connection  with  the  heat  equation: 
starting  with  bounded  initial  data,  f(x),  the  solution  of  the  heat  equation  (2.2.1) 

(3.4.1)  u(x,t)  =  Q*f(x),  Q(x)=-^=e'^r,  t>  0 

V  47T (X 

represents  the  effect  of  smoothing  f(x ),  so  that  u(-,t  >  0 )eC°°  (in  fact  analytic)  and  u(x,t  J. 
0)  =  f(x). 

A  general  process  of  smoothing  can  be  accomplished  by  convolution  with  appropriate 
smoothing  kernel  Q(x) 

(3.4.2)  fe(x )  =  Qe(x )  *  f(x) 
such  that 

(3.4.3a)  Qc(x)*f(x) 

is  sufficiently  smoother  than  f(x)  is,  and 


(3.4.36)  Qe(x)  *  f(x)  — *  f(x). 

e~*0 

With  the  heat  kernel,  the  role  of  e  was  played  by  time  t  >  0.  A  standard  way  to  construct 
such  smoothers  is  the  following.  We  start  with  a  C5-function  supported  on,  say,  (-1,1),  such 
that  it  has  a  unit  mass 


(3.4.4a) 


and  zero  first  r  moments 


(3.4.46)  J  x](f>(x)dx  =  0,  ;  =  l,2,---,r. 

Then  we  set  Qc(x)  =  and  consider 


(3.4.5) 


fc(x)  =  Qe(x)  *  f(x),  e  >  0. 


Now,  assume  /  is  (r  +  1)  -  differentiable  in  the  e  neighborhood  of  x;  then,  since  Qe(x )  is 
supported  on  (— e,e)  and  satisfies  (3.3.4a)  as  well,  we  have  by  Taylor  expansion 


(3.4.6) 


f(x)  -  Qe(x )  *  f(x)  =  f  Qc(y)[f{x)  -  f{x  -  y))dy  = 


(-y)r+1 

(r  -f  1)! 


fr+1(0 


dy. 
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The  first  r  moments  of  Qc{y)  vanish  and  we  are  left  with 

(3.4.7)  \f(x)  -  Qe(x)  *  f(x) |  <  Const.  Max|y_I|<£|/(,"fl)(2/)  •  er+1, 

i.e.,  /£(x)  converges  to  f(x)  with  order  r  +  1  as  e  — >  0.  Morevoer,  /£(z)  is  as  smooth  as  <f>(x ) 
is,  since 

(3.4.8)  /.(*)  =  y'  \q  (^)  f(s)iy 

has  many  bounded  derivatives  as  Q  has,  i.e.,  starting  with  differentiable  function  /  of  order 
r  +  1  in  the  neighborhood  of  x,  we  end  up  with  regularized  function  f:(x)  in  Cs,s  >  r. 


Example:  For  C°°  regularization  -  choose  a  unit  mass  C°°  kernel,  see  Figure  6, 


Figure  6: 


(3.4.96)  |  f(x)  ~  fc(x) |  <  Const. Max|y_x)<£|/,(T/)|  •  e  ->  0. 

To  increase  the  order  of  convergence,  we  require  more  vanishing  moments  which  yield  more 
oscillatory  kernels  as  in  Figure  7. 
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Figure  7: 

We  note  that  this  smoothing  process  is  purely  local  -  it  involves  ^-neighboring  values  of 
GT+1  function  /,  in  order  to  yield  a  ^-regularized  function  fe{x)  with  /c( x) — >-/(x).  The 
convergence  rate  here  is  r  -f  1 . 

We  can  also  achieve  local  regularization  with  spectral  convergence.  To  this  end  set 
(3.4.10a)  Qn(x)  =  p(x)Dx(x) 

where  p(x)  is  a  ^“-function  supported  on  (—£,£)  such  that,  see  Figure  8, 

(3.4.106)  p(  0)  =  1. 


Figure  8: 
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Consider  now: 


(3.4.11) 


/at M  =  Qn  *  /(x)- 


Then  /n{x)  is  because  Qh  is;  and  the  convergence  rate  is  spectral,  since  by  (1.1.25) 


(3.4.12a) 


/(x)  -Qn*  /(x)  =  /(x)  -  f  DN(y)p(y)f(x  -  y)dy 

J\y\<e 

=  /(x)  -  p(y)f(x  -  y ) | y—  o  +  residual 


where 

(3.4.125)  |residual|  <  Const. \\p(y)f(x  -  y)||H*(-e,e)-^7rr  for  any  s  >  0. 


Thus,  the  convergence  rate  is  as  fast  as  the  local  smoothness  of  /  permits;  of  course  with 
p  =  1  we  obtain  the  C^-regularization  due  to  the  spectral  projection.  The  role  of  p  was  to 
localize  this  process  of  spectral  smoothing. 

We  can  as  easily  implement  such  smoothing  in  the  Fourier  space:  For  example,  with  the 
heat  kernel  we  have 


(3.4.13) 


u(k,t)  =  e-akHf(k) 


so  that  u(k,t)  for  any  t  >  0  decay  faster  than  exponential  and  hence  u(x,t  >  0)  belong 
to  Hs  for  any  s  (and  by  Sobolev  imbeddings,  therefore,  is  in  C°°  and  in  fact  analytic).  In 
general  we  apply, 

CO 

(3.4.14)  /.(*)=  £  4(i)/(*=yfe 

k—~oo 

such  that  for  /e(x)  to  be  in  Hs  we  require 

OG 

(3.4.15a)  £  (1  +  \k\2Y\Qc{k)\2{f{k)\2  <  Const. 

k—  —  ct 

and  r  +  1  order  of  convergence  follows  with 

(3.4.156)  \4>e(k)  —  ll  <  Const.(£^)r+1. 


Indeed,  (3.4.15b)  implies 

CO 

I/to  - /.Ml  <  Const.er+1  £  |ir+7(i)e'te| 
(3.4.16)  *=~a 

<  Const.  Max|/r+1)|-£r+1. 
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H- 


Note:  Since  4>c(k)  J.  0  we  can  deal  with  any  unbounded  /  by  splitting  2S|f.|<p.0|  +  S|fc|>|fc0|-  To 
obtain  spectral  accuracy  we  may  use 

[Si,  1*1  <  f 

(3.4.17a)  QN(k)  = 

~  smoothly  decay  to  zero  y  <  |&|  <  N 
Clearly  Q ^  *  f{x)  is  C00  and  as  in  Section  1 

(3.4.174)  \!(x)  -  Qn  .  /(*)!  <  Z  l/(*y“l  <  Const.il/IU.  ■  -A-. 

|A=|>Ar£  iV 

We  emphasize  that  this  kind  of  smoothing  in  the  Fourier  space  need  not  be  local;  rather 
Qe(x)  or  ^at(x)  are  negligibly  small  away  from  a  small  interval  centered  around  the  origin 
depending  on  £  or  jj.  (This  is  due  to  the  uncertainty  principle.) 

The  smoothed  version  of  the  pseudospectral  approximation  of  (3.2.16)  reads 

(3.4.18)  =  TpN{a(x,  t)^{Q  *  v.v)) 

i.e.,  in  each  step  we  smooth  the  solution  either  in  the  real  space  (convolution)  or  in  the  Fourier 
space  (cutting  high  modes).5  We  claim  that  this  smoothed  version  is  stable  hence  convergent 
under  very  mild  assumptions  on  the  smoothing  kernel  Qn(x).  Specifically,  (3.4.18)  amounts 
in  the  Fourier  space,  compare  (3.2.3) 

dv 

(3.4.19)  j-t  =  AnAQnv. 

The  real  part  of  the  matrix  in  question  is  given  by 

(3.4.20a)  [RcAkAQm]^  =  »( A*  -  A ,)  £  a[k  -  j  +  p(2N  +  1)],  -N  <  k,j  <  N 

p 

where  A  Qn  =  diagfc(iA*) 

(3.4.206)  i\k  =  ikQu(k) 

is  interpreted  as  the  smoothed  differentiation  operator.  Now,  looking  at  (3.4.20a)  we  note: 


1.  For  p  —  0  we  are  back  at  the  spectral  analysis,  compare  (3.1.12),  (3.1.13)  and  the  real 
part  of  the  matrix  in  (3.4.20a)  -  the  aliasing  free  one  -  is  bounded. 

2.  We  are  left  with  |p|  =  1:  in  the  unsmoothed  version,  these  terms  were  unbounded 
since  |Afc  —  Ajj  f  co  as  k  j  —  N  or  j  f  N.  With  the  smoothed  version,  these  terms  are 
bounded  (and  stability  follows),  provided  we  have 

5Either  one  can  be  carried  out  efficiently  by  the  FFT. 
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(3.4.21) 


|Ajt  =  ikQtf(k) |  *  0 

Ifcltw 


For  example,  consider  the  smoothing  kernel  Qn(x)  where 

sin  kh  ,  2tt 


(3.4.22a) 


QN(k)  = 


kh 


h  = 


2N  +  1 ' 


This  yields  the  smoothed  differentiation  symbols 


(3.4.226) 


As.  =  fsin 


kh 

T 


which  corresponds  to  the  second  order  center  differencing  in  (1.2.36);  stability  is  immediate 
by  (3.4.21)  for 


_•  2ttA; 
i-i  _  Sm  2AM-1 
I'**  ~  2-irfc 
2W+1 


— >  0 
IMTW 


Yet,  this  kind  of  smoothing  reduces  the  overall  spectral  accuracy  to  a  second  one;  a  fourth 
order  smoothing  will  be 


(3.4.23a) 


or 


(3.4.236) 


At  =  i- 


Afc  = 


kh  .  2 kh 

4  sin  - sm  — — 

h  2  h 


6i  sin  kh 


6,  4  +  2  cos  kh ’  ® N ^  ik' 


In  general,  the  accuracy  is  determined  by  the  low  modes  while  stability  has  to  do  with 
high  ones.  To  entertain  spectral  accuracy  we  may  consider  smoothing  kernels  other  than 
trigonometric  polynomials  (=  finite  difference),  but  rather,  compare  (3.4.17) 

=  1,  |*|<f 


(3.4.24) 


QN(k)  = 


^  smoothly  decay  to  zero  ~  <  |&|  <  N. 

An  increasing  portion  of  the  spectrum  is  differentiated  exactly  which  yields  spectral  accuracy; 
the  highest  modes  are  not  amplified  because  of  the  smoothing  effect  in  this  part  of  the 
spectrum. 

We  close  this  section  noting  that  if  some  dissipation  is  present  in  the  differential  model 
to  begin  with,  e.g.,  with  the  parabolic  equation 


(3.4.25) 


du  _  d 
dt  dx 


a(x ,  t) 


du N 
dx  , 


a(x,  t)  >  a  >  0, 


then  stability  follows  with  no  extra  smoothing.  The  parabolic  dissipation  compensates  for 
the  loss  of  “one  derivative”  if  first  order  terms  are  present.  To  see  this  we  proceed  as  follows: 
multiply 


(3.4.26) 


dyN ,  ,x  d 

-dTM=K 


t)  0 
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by  v^(xut)  and  sum  to  obtain 


(3.4.27) 


2  dt  ?  Q  ~  S  M*v,  i) (**,  *))• 


Suppressing  excessive  indices,  we  have  for  the  RHS  of  (3.4.27) 

?• 4  («*>£)-  \  e  £  {■■(»’£) -  ? (£)' 1 . 

Now,  the  first  sum  on  the  right  gives  us  the  usual  loss  of  one  derivative  and  the  second 
are  compensates  with  gain  of  such  quantity.  Petrowski  type  stability  (gain  of  derivatives) 
follows.  We  shall  only  sketch  the  details  here.  Starting  with  the  first  term  on  the  right  of 


(3.4.28)  we  have 


(3-4-29)  2  ^  dx  {avVv  dx  )  ~  2  /  ‘1  +  taliasin§  errors] 


while  for  the  second  term 


J2av(t) 


(3-430)  -£•*>(£)  W  [£<*•* )f* 

and  this  last  term  dominates  the  RHS  of  (3.4.29). 
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APPENDIX 


A.l.  Fourier  Collocation  with  Even  Number  of  Gridpoints 


We  assume  w(x)  is  known  at 


(Al.l) 


=  w(xv)  x„  =  r  4-  vh  v  —  0, 1,  •  •  •  f2N 


with  h,  =  =  J?  and  0  <  r  <  h.  We  use  the  trapezoidal  nodes 


(A.1.2) 


1  2  N  1  2N-1 

m  =  7-  E  'V“"‘  =  777  E 


to  obtain  the  pseudospectral  approximation 


(A.l. 3) 


i>Nw=  £  "w{k)eikx. 


Note:  We  now  have  only  2 N  pieces  of  discrete  data  at  the  different  2 N  grid  points 
x0,  Xj,  •  •  • ,  Z2W-1  and  they  correspond  to  2 N  waves,  as  we  have  a  “silent”  last  mode,  i.e., 
with  r  =  0,  k  =  N,  Im[e'kx }X-Xl/  =  •£  sin  =  0.  This  is  a  projection,  since  in  view  of  (A. 1.3) 
ipNW  is  the  interpolant  of  w(x)  at  x  =  xv\ 


(A.l. 4) 


N  i  2W+1  ] 

=  E  "  777  E 

k=-N  lAiV  V=0 


2N-\  i  N 

=  E  "’(*-)  •  777  E  =  »(*,.)• 

u=0  ZJV  k=-N 


The  aliasing  relation  in  this  case  reads  -  compare  (1.2.7) 


(A.1.5) 


w(k)=  Y1  e,p2Nrw(k  +  2  PN) 


and  spectral  convergence  follow  -  compare  with  (1.2.16) 


(A. 1.6) 


||A^to(x)|!tf.  <  Const,  •  ||rAru;(x)|jH-)  5  >  -. 


In  the  usual  sin-cos  formulation  it  takes  the  form 
(A.  1.7) 

N  _  r  h  1  2N+1  t, 

"ak c°s kx  +  h sin kx,  ~k  =  —  ^  w(xv)  cos  Xv  }  0<k<N. 

k= o  °k  N  u=0  sin  kxv 

Noting  that  =  0  we  have  2 N  free  parameters  a0,  {a*,  ^>k}k=i  and  a ^  to  match  our  data 
at 
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